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1.0  INTRODUCTION 


This  report  contains  a  summary  of  progress  on  AFOSR  contract  F49620-79-C-0057 
on  nonlinear  finite  element  analysis  and  outlines  suggested  further  research. 
Section  2.1  begins  by  describing  the  deficiencies  in  the  state-of-the-a~t  in 
nonlinear  analysis  as  they  relate  to  and  have  motivated  the  present  research. 
Section  2.2  discusses  methodology  and  progress  in  the  current  contract  in 
detail.  Section  2.3  presents  results  for  several  large  deflection  problems 
for  beams,  addressing  features  of  the  nonlinear  behavior  which  illustrate 
both  the  advantages  and  the  deficiencies  of  the  developed  methods.  Section 
2.4  suggests  areas  for  future  studies  related  to  the  work.  Appendices  A  and  B 
provide  technical  details  on  two  areas  of  development  which  are  briefly  dis¬ 
cussed  in  Sections  2.1  and  2.2,  but  for  which  the  reader  may  desire 

« 

clarification.  Appendix  C  reproduces  portions  of  the  proposal  (Reference  5) 
for  the  present  contract,  for  convenient  reference  in  the  present  context. 

The  goal  of  the  present  research  has  been  the  development  and  evaluation  of 
improved  displacement-method  finite  element  approaches  for  the  analysis  of 
structural  problems  with  geometrical  nonlinearities.  The  initial  work  In 
this  subject  was  done  by  Haftka,  Mallett  and  Nachbar  in  reference  1.  Further 
work  by  Jones  (References  2,  3,  4)  followed  along  the  lines  outlined  in 
reference  1.  Reference  3  concluded  that  an  improved  approach  could  be  devel¬ 
oped  through  new  finite  element  formulations  coupled  with  a  stepwise  non¬ 
linear  solution  procedure.  Reference  5  proposed  and  outlined  the  development 
of  these  new  approaches,  which  have  been  pursued  In  the  present  AFOSR- 
sponsored  research  and  are  reported  herein. 

There  are  three  technical  areas  associated  with  geometrically  nonlinear 
finite  element  analysis  which  require  investigation  in  order  to  achieve  the 
desired  research  results.  These  are  introduced  briefly  here  in  order  to  put 
into  perspective  more  detailed  discussions  in  the  sections  which  follow. 

First,  it  is  required  to  derive  a  new  type  of  finite  element  which  is 
numerically  well  behaved  when  the  total  nonlinear  deformations  are  retained 
in  analysis.  This  requirement  addresses  the  role  of  the  nonlinear  contribu¬ 
tions  to  the  stresses  as  major  controlling  factors  in  the  equilibrium  state. 
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With  conventional  elements,  retaining  these  stresses  results  In  serious 
errors  with  consequent  Incorrectly  convergent,  or,  frequently,  divergent 
solution  calculations.  Secondly,  It-ls  required  to  develop  geometrical 
representations  of  the  deformations  which  avoid  any  cumulative  Inconsistency 
between  the  deformations  and  the  displacements.  Such  errors  cannot  be  elimi¬ 
nated  by  Iteration.  Finally,  it  is  required  to  develop  improved  stepwise 
solution  procedures,  with  residual  load  iteration,  which  are  convergent  for 
large  load  steps.  The  presence  of  large  residual  loads  caused  by  the 
nonlinear  deformations  often  causes  failure  of  conventional  solution 
procedures.  These  three  requirements  address  proven  deficiencies  in  the 
application  of  conventional  finite  element  metnods  to  strongly  nonlinear 
analysis,  as  manifested  by  both  unacceptable  Inaccuracies  In  problem  solu¬ 
tions  and  serious  difficulties  In  obtaining  convergent  stepwise  solution 
processes. 


2.0  TECHNICAL  DISCUSSION 


This  section  discusses  the  goals  of  the  present  research,  progress  on  the 
current  contract,  and  recomnended  areas  for  further  research. 

2.1  Goals  and  Approach  of  Current  Research 

The  physical  problems  which  have  motivated  the  current  research  primarily 
Involve  those  structures  In  which  the  nonlinearly  induced  stresses  due  to 
large  rotations  are  critical  in  determining  the  correct  equilibrium  state  cf 
the  structure.  Examples  of  such  problems  include  such  difficult  problems  as 
the  buckling  of  shells  and  the  postbuckllng  action  of  panels,  and  also  simple 
cases  such  as  the  stretching  of  cables  and  bending  of  beams.  In  all  of  these 
case?  it  is  necessary  to  include  a  complete  and  accurate  description  of  the 
nonlinear  stress  field  within  the  structure  in  order  to  correctly  address  the 
equilibrium  problem.  In  general,  finite  element  analysis  has  encountered  a 
great  deal  of  difficulty  in  doing  this.  The  difficulties  arise  because 
conventional  element  formulations  retaining  complete  nonlinear  stress  calcu¬ 
lations  encounter  several  fundamental  problems.  One  of  these  problems,  and 
that  which  requires  the  development  of  new  types  of  elements,  is  that  the 
nonlinearly  induced  stresses  include  physically  unrealistic  components  which 
become  "locked"  within  the  elements,  unremovable  except  by  major  reduction  of 
the  displacement  magnitudes.  This  causes  excessive  structural  stiffness  and 
results  in  inaccurate  problem  solutions.  Another  problem  is  that,  due  to 
including  the  nonlinear  stresses,  the  solution  procedure  is  required  to  deal 
with  very  large  residual  loads  in  the  iterative  portion  of  the  calculation 
process.  It  is  usually  found  that  such  large  residual  loads  cause  solution 
divergence  or  excessively  slow  convergence  and  very  large  computational 
costs.  This  difficulty  requires  the  development  of  solution  procedures  which 
are  improved  in  character  and  tailored  specifically  to  problems  having  large 
nonlinear  stresses. 

An  additional  difficulty  is  that  in  problems  of  the  types  under  discussion 
there  are  generally  large  rotations  and  corresponding  deformations  which  are 
computed  in  a  stepwise  manner  during  the  solution  process.  Total  deformations 


are  often  determined  by  sunmlng  increments.  This  approach  can  cause  cumula- 
tive  errors  in  the  solution  which  are  not  correctable  through  the  residual 
load  iteration  process.  There  are  several  sources  of  such  errors.  The  first 
and  that  most  commonly  encountered  In  conventional  approaches  is  a  cumulative 
error  in  the  strains  caused  by  the  determination  of  the  strain  through  step- 
wise  Incrementation.  This  error  results  from  linearization  of  the  strain 
Increments,  and  Is  usually  quite  large.  A  second  source  of  error  occurs  In 
cases  having  moderately  large  rotations  in  three  dimensions.  The  representa¬ 
tion  of  such  three-dimensional  rotations  through  incrementation  of  cartesian 
rotations  will  cause  both  Incorrect  total  nonlinear  strains  and  also  a  cumula¬ 
tive  error  In  the  orientation  of  the  structure  in  space.  The  error  will  make 
the  strains  inconsistent  with  the  true  total  displacement  and  rotation  state, 
and  hence  is  not  recoverable  through  residual  load  iteration.  This  rotation 
problem  must  be  solved  In  order  to  develop  finite  element  procedures  which  are 
applicable  to  such  problems  as  the  combined  lateral  and  torsional  nonlinear 
deformation  of  beam  structures.  Another  type  of  error  occurs  in  total  Lagran- 
glan  formulations,  and  results  from  an  exchange  of  roles  between  the  bending 
and  membrane  displacements  when  the  rotations  become  large. 

The  above  discussions  point  clearly  to  three  primary  technical  goals  for  the 
current  research.  The  first  Is  the  development  of  new  element  types  which  are 
formulated  sped fl cal ly  for  including  complete  nonlinear  strains  In  finite 
element  calculations.  The  beam  and  shell  elements  under  development  in  this 
research  accomplish  this  through  the  use  of  interior  modes  to  Incorporate 
higher  order  axial  (beam)  and  membrane  (plate,  shell)  displacement  functions. 

The  second  primary  goal  Is  to  avoid  unrecoverable  cumulative  error,  that  Is, 
any  type  of  error  which  Is  unrecognized,  and  hence  uncorrectable,  by  the 
residual  load  Iteration  process.  The  cumulative  errors  due  to  incrementing 
element  strains  In  a  stepwise  procedure  have  beer  avoided  by  computing  the 
total  nonlinear  strains  directly,  rather  than  by  Incrementation.  A  total 
Lagrangian  formulation  with  updating  is  used  to  accomplish  this.  The  cumula¬ 
tive  error  In  the  structural  rotation  state  due  to  summing  cartesian  rotation 
Increments  has  beer  avoided  through  a  rigorous  three-dimensional  rotation 
description.  In  this  approach,  the  rotation  state,  both  total  and 
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Incremental,  Is  represented  by  three  sequenced  angles,  called  Euler  angles, 
which  uniquely  define  the  total  rotation  state  for  arbitrarily  large  mutlcns. 
This  appears  to  be  a  new  approach  In  finite  element  formulations. 

The  third  primary  goal  Is  technically  distinct  from  the  finite  element 
research  of  the  work,  but  was  necessary  In  order  to  achieve  nunerlcal  verifi¬ 
cations  of  the  element  technology.  It  Is  the  development  of  solution  proced¬ 
ures  which  are  rapidly  convergent  despite  the  nunerlcal  difficulties  Inherent 
In  the  solution  process  for  strongly  nonlinear  problems.  Several  approaches 
have  been  Investigated  In  this  regard.  The  Initial  approach  used  an 
Internally  nonlinear  stepwise  solution  procedure  with  Iteration  of  the  resid¬ 
ual  load  state  after  each  load  step.  The  nonlinear  stepwise  capability  was 
based  on  the  static  perturbation  procedure  (References  6,  7,  8).  Additional 
means  to  convergence  acceleration  were  found  necessary  In  conjunction  with 
the  static  perturbation  method;  these  are  discussed  in  Section  2.2.  Ulti¬ 
mately,  this  approach  was  found  Inadequate  (except  for  small  load  steps  and 
moderately  small  displacements),  and  a  method  utilizing  adaptive  modification 
of  the  structural  stiffness  matrix  (the  BFGS  method,  reference  10)  was  suc¬ 
cessfully  implemented  in  place  of  the  static  perturbation  method. 

The  total  set  of  technical  goals  for  the  present  research  include.  In  addition 
to  the  three  primary  goals  described  above,  a  number  of  important  secondary 
Items.  These  have  been  grouped  into  two  categories.  The  first  category 
relates  to  the  matter  of  suitable  finite  element  strain  displacement  equa¬ 
tions.  It  is  required  that  the  basic  rules  governing  finite  element  displace¬ 
ment  states  and  strain-displacement  relations  be  followed:  the  element  can 
undergo  rigid  body  displacements  and  constant  strain  states;  large  rigid  body 
motions  must  not  cause  element  deformations.  These  requirements  are  met  by 
using  cartesian-based  elemental  coordinate  systems.  Other  requirements 
relate  to  the  degree  of  approximation  of  the  nonlinear  portions  of  the  strain- 
displacement  equations  and  the  accurate  representation  of  geometry  for  shell 
elements.  It  Is  necessary  to  keep  the  strain-displacement  equations  simple 
and  amenable  to  nunerlcal  processing  as  large  displacements  and  rotations 
develop.  In  particular.  It  Is  required  to  avoid  the  complexity  of  a  nonlinear 
shell  strain  displacement  formulation  of  the  intrinsic  coordinate  type.  The 
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use  of  the  Mar  guerre  type  of  strain  formulation,  together  with  element  local 
coordinate  system  updating,  has  fulfilled  these  requl ranents .  The  second 
category  relates  to  generality  of  solution  procedures.  It  was  originally 
Intended  that  the  solution  procedures  developed  In  this  work  be  applicable  to 
an  extended  set  of  of  problem  types  Including  buckling  (bifurcation,  limit 
points  and  snap  through),  and  also  that  the  solution  procedure  be  extendable 
to  the  case  of  dynamic  response  calculations  where  strong  nonl Inearl  ties  are 
present.  The  static  perturbation  approach  appeared  to  have  the  potential  to 
satisfy  these  requirements.  The  difficulties  in  achieving  convergence  of 
nonlinear  stepwise  solutions  with  the  static  perturbation  method,  and  the 
adoption  of  the  3FGS  method  for  this  purpose,  have  necessarily  modified  the 
original  plan  to  develop  a  solution  procedure  with  direct  applicability  to 
both  dynamics  and  buckling  problems. 

2.2  Progress  of  Current  Research 

In  the  work  accomplished  to  date  the  theoretical  development  of  the  finite 
element  formulations  for  the  new  type  of  two-dimensional  and  three-dimens¬ 
ional  beam  elements  has  been  completed  and  the  Euler  angle  theory  for  both  the 
beam  elements  and  the  plate  and  shell  elements  has  been  developed.  By  calcu¬ 
lations,  It  has  been  demonstrated  that  the  type  of  element  formulation  under 
development  Is  completely  successful  In  handling  the  large  stresses  inherent 
in  the  large  rotation  nonlinear  state.  A  technically  advanced  set  of  solution 
procedure  algorithms  has  been  developed,  refined  and  verified  through  numeri¬ 
cal  studies  using  a  two-dimensional  beam  element  code.  The  solution  proced¬ 
ures  appear  to  be  rapidly  convergent  for  large  step  sizes  and  strong  nonl Ine¬ 
arl  ties.  A  number  of  technical  difficulties,  seme  expected  and  some  unantici¬ 
pated,  have  been  encountered.  The  discussion  which  follows  attempts  to 
describe  the  chronology  of  the  work  accomplished  and  its  present  status, 
covering  each  task  and  difficulty  In  some  detail. 

The  current  research  began  with  the  development  of  the  element  technology  for 
a  curved  beam  element  whose  Initial  shape  and  subsequent  deformation  are 
constrained  to  take  place  in  a  single  plane.  This  has  been  called  the  2-0 
beam  element  (see  Appendix  C,  Section  2.2  and  Figure  1).  The  element  has  five 
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nodes:  nodes  1  and  5  are  the  end  nodes;  2,  3,  and  4  are  internal.  In  its 
present  form,  the  nodes  are  equally  spaced  along  the  length  of  the  element. 

Nodes  1,  3  and  5  each  have  three  freedoms.  Including  the  axial  displacement, 
the  lateral  or  bending  displacement,  and  the  bending  rotation.  Nodes  2  and  4 
have  only  the  axial  displacement  freedom.  This  unique  nodal  and  freedom 
arrangement  provides  axial  displacements  which  are  quartic  functions  and 
bending  displacements  and  rotations  which  are  quadratic.  By  this  means  accur¬ 
ate  stress  and  strain  representations  are  obtained  despite  the  rotations 
present  in  strongly  nonlinear  problems  (see  discussion  in  Appendix  C,  Section 
2.2,  Stability  Elements).  The  development  of  this  element  encountered  sev¬ 
eral  technical  difficulties  related  to  its  unusual  nodal  and  freedom  arrange¬ 
ment  and  in  regard  tr  geometrical  updating.  For  example,  In  performing  a 
rotational  updating  transformation,  the  calculation  of  the  transformed  values 
of  the  axial  freedoms  at  nodes  2  or  4  requires  including  in  the  transformation 
the  potentially  large  bending  displacements  at  these  nodes.  Since  the  bending 
displacement  Is  not  an  available  freedom  at  rodes  2  and  4,  Its  value  must  be 
generated  by  interpolation  using  the  bending  displacement  values  at  nodes  1,  3 

and  5.  Related  difficulties  are  encountered  in  transforming  the  loading  on 

% 

the  beam  element.  The  coupling  together  of  the  axial  and  bending  displace¬ 
ments  (or  loads)  in  such  updating  transformations  Is  the  means  by  which  any 
small  angle  approximations  used  In  the  nonlinear  strain  equations  are  removed 
by  updating.  This  has  implications  regarding  the  performance  of  solution  proce¬ 
dures,  the  calculation  of  residual  loads,  and  the  criteria  controlling  updating. 

The  next  step  in  the  research  was  the  development  of  solution  procedures 
appropriate  to  the  implementation  of  the  2-D  beam  element.  The  solution 
procedure  development  was  based  on  the  static  perturbation  method  of  the 
second  order  (the  static  perturbation  method  is  discussed  in  Appendix  C, 

Section  2.2,  Nonlinear  Step  Static  Solution  Procedure,  and  also  in  a  mathe¬ 
matical  derivation  starting  on  page  C30  of  Appendix  C.  See  also  Appendix  A, 
page  AT,  for  a  discussion  of  path  parameters).  In  this  formulation,  the 
structural  displacements  are  expressed  as  a  second  degree  Taylor  series  in  a 
path  parameter.  The  path  parameter  Is  the  Taylor  series  argument.  In  the 
initial  development  of  this  theory  the  path  parameter  was  the  load  itself. 
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That  is,  the  structural  displacements  were  expressed  as  a  second  order  Taylor 
series  in  the  load.  This  approach  Is  computationally  simple  and  provides 
excellent  results  for  many  types  of  problems.  It  presumes  that  the  displace¬ 
ments  are  well  behaved  functions  of  the  load,  and  in  particular  that  displace¬ 
ments  cannot  occur  without  corresponding  changes  in  the  value  of  the  load. 
Hence  this  approach  would  be  unsuitable  for  buckling  problems,  in  which  dis¬ 
placements  can  occur  at  constant  or  nearly  constant  load. 

A  proof-of-concept  computer  program  was  written  to  Implement  this  procedure 
In  conjunction  with  the  2-D  beam  element.  The  program  was  designed  to  have 
maximum  adaptability  to  future  extensions  of  element  technology.  In  numeri¬ 
cal  work  with  this  computer  program,  the  surprising  result  was  found  that  for 
certain  problem  types  (large  rotations  with  very  small  axial  stresses)  the 
second  order  static  perturbation  procedure  often  displayed  poor  convergence 
or  even  divergence  unless  the  applied  load  steps  were  made  small.  In  an 
attempt  to  improve  this  situation  the  static  perturbation  procedure  was 
extended  to  include  the  third  degree  Taylor  series  terms.  This  extension  led 
to  a  great  deal  of  complexity,  both  in  the  area  of  theory  development  and  also 
In  coding  work,  and  produced  disappointing  results.  In  nunerical  work  it  was 
found  that  the  second  order  approach  consistently  performed  better  than  the 
substantl ally  more  complex  third  order  formulation.  A  considerable  amount  cf 
study  was  done  to  explain  this  unexpected  result  and  to  understand  more  fully 
the  behavior  of  the  second  and  third  order  formulations. 

It  was  determined  that  the  second  order  procedure  provides  corrective  dis¬ 
placements  which  primarily  reduce  errors  In  the  axial  force  equilibrium 
state.  The  third  order  process,  on  the  other  hand,  primarily  provides  correc¬ 
tive  bending  displacements  In  order  to  reduce  errors  in  the  bending  (lateral) 
load  equilibrium  state.  The  static  perturbation  approaches  accomplish  these 
corrections  through  a  type  of  residual  load  evaluation,  which  is  made  using 
only  the  start-of-step  geometry  and  deformation  description  of  the  structure. 
Since  this  start-of-step  state  is  approximate  in  its  ability  to  forecast  end- 
of-step  residual  loads,  a  corresponding  approximation  occurs  in  the  static 
perturbation  corrective  displacement  values.  In  contrast,  corrective 
displacements  computed  by  the  stepwise  Iterative  process,  using  rigorous. 
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end-of-step  evaluated  residual  loads,  are  able  to  respond  exactly  to  the 
nonlinearity- induced  errors  which  have  occurred  during  a  given  load  incre¬ 
ment.  In  geometrically  nonlinear  analysis,  axial  equilibrium  corrections  are 
crucially  Important,  because  the  axial  forces  combine  with  the  rotations  to 
produce  potentially  large  bending  equilibrium  errors.  Hence  the  second  order 
static  perturbation,  despite  its  approximate  nature,  is  beneficial  to  solu¬ 
tion  convergence.  On  the  other  hand,  errors  in  bending  displacement  predic¬ 
tion  can  be  crucially  damaging  to  solution  convergence,  because  of  their 
potential  to  cause  large  rotations,  nonlinear  axial  strains,  and  hence  large 
axial  force  errors.  Hence  the  approximations  in  the  bending  displacement 
corrections  computed  by  the  third  order  static  perturbation  process  appear 
unacceptable.  It  was  concluded  for  this  reason  that  conventional  static 
perturbation  of  the  second  order  is  superior  to  the  third  order  procedure  for 
geometrically  nonlinear  analysis  of  "thin"  structures  (beams,  plates, 
shells) . 

At  this  point  in  the  research  it  was  felt  worthwhile  to  extend  the  static 
perturbation  approach  to  a  more  general  type  of  formulation,  in  which  the 
Taylor  series  path  parameter  is  def omation-rel ated  (Appendix  A  discusses  the 
path  parameter  in  detail).  A  useful  path  parameter  of  this  type  is  similar  to 
the  structural  strain  energy  function,  taking  the  form 

S2  =  &CT  K 

where  S  is  the  path  parameter,  K  is  the  tangent  stiffness  matrix,  and  iQ  is 
the  incremental  displacement  of  any  load  (cr  iteration)  step.  This  type  of 
path  parameter  has  the  advantage  of  applicability  to  buckling  problems.  It 
was  felt  that  calculations  using  this  particular  type  of  path  parameter  might 
shed  some  light  on  the  overall  behavior  of  the  static  perturbation  process  In 
the  types  of  problems  under  study.  This  particular  extension  again  led  to  a 
great  deal  of  complexity,  in  both  theory  and  numerical  approach.  Nunerlcal 
work  with  this  approach  was  again  disappointing,  and  showed  that  for  problems 
in  which  the  displacements  are  well  behaved  functions  of  the  load,  the  energy- 
based  path  parameter  formulation  does  not  provide  any  advantages  over  the 
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load-based  path  parameter  approach.  Only  in  cases  in  which  the  load-displace¬ 
ment  relationship  is  poorly  behaved,  as  occurs  in  buckling  problems,  would  the 
generality  of  this  approach  be  advantageous. 

In  order  to  obtain  Improved  solution  procedure  performance,  recourse  was  made 
to  methods  developed  during  previous  experience  (prior  AFOSR  contract,  refer¬ 
ences  3,  4)  in  nonlinear  analysis  of  beams  and  plates.  In  this  work  it  had 
been  concluded  that  the  axial  force  equilibrium  errors  are  primarily  a  result 
of  bending  displacements  which  take  place  without  perfectly  "matched"  axial 
displacements.  The  axial  force  errors  are  usually  very  large.  Together  with 
the  rotations  they  cause  large  error  loads  acting  on  the  relatively  flexible 
bending  displacements.  This,  in  turn,  causes  further  bending  displacement 
errors  with  even  larger  axial  force  errors.  Thus,  the  errors  in  the  axial 
forces  have  a  tendency  to  magnify  themselves  and  cause  divergent  calcula¬ 
tions.  In  the  previous  work  it  was  found  that  a  successful  method  of  accele¬ 
rating  convergence  is  to  perform  residual  load  iterations  for  the  single 
purpose  of  "matching"  the  axial  displacements  to  the  bending  displacements, 
thus  removing  the  axial  load  equilibrium  errors.  This  is  implemented  by  a 
solution  procedure  employing  "alternate-freedom"  iterative  corrections.  In 
this  procedure,  the  first  incremental  displacement  of  a  load  step  includes  all 
of  the  freedoms  of  the  finite  element  model.  The  next  increment  Is  the  first 
iteration.  It  only  includes  the  axi  al/ membrane  freedoms,  and  thereby  relaxes 
the  axial/membrane  force  errors.  The  next  Increment  is  the  second  iteration; 
it  Includes  all  freedoms.  The  third  iteration  Includes  only  the  axial/mem¬ 
brane  freedoms,  and  so  on.  It  has  been  found  that,  despite  element  curvature 
and  prior  deformation,  in  each  axi  al /membrane- freedom- only  Iterative  correc¬ 
tion  the  axial/membrane  force  equilibrium  errors  are  significantly  reduced. 
This  prevents  these  errors  from  causing,  in  the  subsequent  iteration,  large 
bending  displacement  errors.  This  type  of  solution  procedure  always  accele¬ 
rates  convergence,  and  often  achieves  convergent  solutions  where  other 
approaches  encounter  divergence.  The  alternate-freedom-iteration  procedure 
was  added  to  the  static  perturbation  method  to  achieve  a  combined  process 
having  the  benefits  of  both  procedures.  The  extension  was  accomplished  such 
that  both  the  all-freedom  and  the  axi al/membrane  freedom  Iterations  are  per¬ 
formed  by  the  static  perturbation  procedure. 
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In  calculations  done  with  this  method,  it  was  found  that  the  bending  moment 
residual  loads  tended  to  remain  excessive  after  the  axial/membrane  and  trans¬ 
verse  load  residuals  were  reduced  to  relatively  small  values.  Consequently, 
the  alternate-freedom-iteration  procedure  was  modified  to  include  both  the 
axial  /membrane  freedoms  and  the  rotational  freedoms.  In  this  form  the  solu¬ 
tion  procedure  appeared  optimum. 

The  above- discussed  solution  procedure  functioned  well  for  problems  with 
moderate  displacements  and  rotations,  but  diverged  when  very  large  displace¬ 
ment;;  (e.g.,  half  the  length  for  a  simple  end-loaded  cantilever)  were  compu¬ 
ted.  To  prevent  the  divergence,  a  procedure  called  a  "line-search"  was 
implemented.  In  this  method,  the  amplitude  of  a  computed  displacement  incre¬ 
ment  Is  scaled,  or  optimized,  in  such  a  way  as  to  minimize  the  solution 
errors,  as  measured  by  residual  load  magnitudes,  which  correspond  to  the  total 
displacements  at  the  end  of  the  Increment.  This  avoids  the  use  of  a  computed 
Increment  where  that  Increment  would  Increase,  rather  than  decrease,  the 
residual  loads.  The  implementation  involves  evaluating  a  measure  of  the 
error,- e.g.,  the  root-sun-square  of  the  residual  loads,  for  several  ampli¬ 
tudes  of  the  computed  Increment,  and  interpolating  on  the  amplitude  to  obtain 
a  minimum  error.  The  Interpolation  Includes  the  axi al/membrane/ rotational 
iteration  corrections.  This  approach  performed  well,  particularly  when 
employed  with  judicious  updating  of  the  structural  stiffness  matrix  in  order 
to  assure  that  the  "shape"  of  the  computed  Increment  Is  a  good  one.  If  the 
stiffness  matrix  Is  not  updated.  It  sometimes  occurs  that  the  relative  magni¬ 
tudes  of  the  incremental  values  of  the  structural  freedoms,  i.e.,  the  "shape" 
of  the  Increment,  Is  sufficiently  Inaccurate  that  even  a  near  zero  amplitude 
of  the  computed  Increment  will  Increase  the  error  level.  In  this  case  the 
line-search  falls,  and  an  accurate  problem  solution  is  not  obtained. 

To  avoid  this  difficulty,  a  method  based  on  stepwise  modification  of  the 
structural  stiffness  mat-'x  can  be  used.  Such  an  approach  involves  stepwise 
modifications  of  the  sJ  -.ness  matrix  such  that  the  matrix  to  be  used  for  the 
next  Increment  Is  the  one  which,  had  It  been  used  for  the  1  ast  Increment, 
would  have  produced  an  accurate  Incremental  nonlinear  response  to  the  step- 
wise  Incremental  loads.  The  procedure  implemented  was  the  BFGS  method, 
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described  in  reference  10.  The  BF6S  method  requires  a  line-search  for  each 
increment;  the  above-described  line-search  method  appears  suitable  and  was 
used,  though  It  differs  somewhat  from  the  one  discussed  in  reference  10.  The 
use  of  the  axi al/membrane/rotational  iterations,  described  earlier,  was  found 
to  be  still  necessary  for  convergence,  and  this  procedure  was  combined  with 
the  BFGS  approach  as  follows.  Each  iterative  increment  was  corrected  by  an 
axial/membrane/rotational  iteration,  and  the  “double  increment"  thus  computed 
was  subjected  to  a  line-search.  The  axi al /membrane/ rotatidnal  correction  was 
separately  computed  for  each  interpolation  amplitude  of  the  line-search,  in 
order  to  account  properly  for  the  effects  of  nonlinearities.  Without  these 
corrections,  the  line-search  always  fails.  The  residual  loads  resulting 
after  the  increment  are  compared  with  those  imposed  at  the  start  of  the 
increment,  resulting  in  the  transformation  matrix  of  the  BFGS  method. 

The  solution  procedure  described  above  was  found  to  be  uniformly  convergent 
for  both  large  and  small  displacements  and  for  large  load  step  sizes.  I;i 
addition,  the  convergence  was  found  to  be  rapid  (typically  15  or  less  itera¬ 
tions).  However,  the  convergence  limit  is  to  a  nonzero  error  level  which 
cannot  be  improved  upon  by  further  Iterations.  This  minimum  error  level  is  a 
function  of  the  amplitude  of  the  displacements  and  also  appears  to  be  influ¬ 
enced  by  the  transformations  associated  with  geometrical  updates  of  the  ele¬ 
ment  baseplanes.  For  the  case  of  a  simple  cantilever  beam,  good  accuracy  is 
obtained  when  the  end  displacement  is  less  than  about  half  the  length  of  the 
beam.  This  appears  to  be  true  whether  one  or  more  elements  are  used  in  the 
finite  element  model.  The  intentions  of  the  current  work  are  to  handle  much 
larger  displacements  than  this  apparent  limitation. 

A  nunber  of  numerical  experiments  were  carried  out  in  attempts  to  understand 
and  eliminate  the  minlmim  error  level  problem.  These  included:  comparisons 
between  results  when  the  transverse  shear  and  extenslonal  deformations  are 
represented  by  strain  equations  permitting  large,  as  opposed  to  moderate, 
rotations;  elimination  (in  deformation  and  residual  loads  calculations)  of  the 
quadratic- In-x  component  of  the  transverse  shear  strain;  the  use  of  double¬ 
precision  arithmetic  In  the  geometrical  updating  calculations;  computing 
large  deflection  solutions  with-and-without  geometrical  and  stiffness  matrix 
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updates.  The  numerical  experiments  involving  the  large  rotation  strain 
equations  improved  but  did  not  eliminate  the  minimum  error  level  problem. 

The  use  of  double  precision  geometrical  update  calculations  resulted  in 
somewhat  modified  solutions,  indicating  numerical  sensitivity  in  this  type 
of  calculation.  This  suggests  the  probable  need  for  extension  to  double 
precision  arithmetic  in  all  stress,  deformation,  and  load  computations. 
However,  the  minimum  error  level  was  not  appreciably  affected  by  this 
experiment. 

Through  omission  of  geometrical  updating,  it  was  found  that  the  minimum  error 
level  could  be  made  essentially  zero.  However,  this  is  not  a  satisfactory 
solution  to  the  problem,  because  it  leads  to  a  non-updated  total  Lagrangian 
approach  which  is  subject  to  important  limitations  and  errors.  In  particular, 
the  desirable  emission  of  certain  normally  negligible  terms  in  the  nonlinear 
strain  equations  is  not  admissible  for  a  non-updated  total  Lagrangian 
approach.  The  geometrical  updating  used  in  the  present  approach  is  considered 
a  valuable  feature,  not  to  be  omitted  as  a  solution  to  the  minimum  error 
problem.  What  has  been  gained  through  the  geometrical  updating  experiment  is: 
it  has  been  verified  that  the  element  itself  is  capable  of  a  "zero"  error 
level  in  the  large  displacement  state;  it  is  clearly  indicated  that  the 
geometrical  updating  introduces  displacement  and  deformation  forms  which  for 
some  reason  are  not  satisfactorily  handled  by  the  iterative  solution  proce¬ 
dure.  It  appears  that  the  transverse  shear  strain  becomes  "locked"  in  the 
element,  causing  the  nonzero  minimum  error  level.  Whether  the  fault  lies  in 
the  element  itself,  or  in  the  solution  procedure,  is  not  clear.  It  appears 
that  the  element  should  provide  convergence  to  a  correct  shear  strain  and 
stress,  since  without  geometrical  updating  it  does  so,  and  also  since  its 
nodal  and  freedom  arrangements  and  its  elemental  coordinate  system  satisfy 
the  required  rigid  motion  and  constant  strain  conditions.  On  the  other  hand, 
recent  literature  suggests  that  similar  elements  (though  without  the  same 
internal  nodes  and  freedoms)  may  have  problems  of  a  similar  nature  to  those 
encountered  in  the  present  work  (References  11,  12).  Thus,  Ine  possibility  of 
an  element-level  problem  cannot  yet  be  dismissed. 


The  total  requirements  for  the  solution  procedure  Include  a  number  of  options 
In  addition  to  those  of  the  basic  Iteration  process  discussed  above.  Other 
needed  features  include:  conditional  geometrical  updating  of  element  local 
coordinate  systems  and  displacement  states;  conditional  updating  of  the 
stiffness  matrix,  with  separate  handling  of  the  linearized  portion  of  the 
matrix  and  the  stress-dependent  ("geometric  stiffness")  portion;  conditional 
controls  on  solution  continuance  or  termination;  limits  on  the  ntmber  of 
geometrical  and  stiffness  updates  and  on  the  error  measure  which  constitutes 
convergence. 

The  developed  solution  procedure  contains  user  options  for  controlling  all  of 
these  items.  Table  1  gives  a  brief  simmary  of  the  total  set  of  solution 
procedure  options.  The  static  perturbation  controls  are  covered  in  part  (a), 
and  the  controls  over  the  entire  procedure  (as  currently  coded)  are  given  in 
part  (b).  One  ’tern  which  requires  further  discussion  Is  the  conditional 
updating  of  the  stiffness  matrix.  The  total  stiffness  matrix  is  the  sun  of 
the  basic  stiffness  matrix  and  the  geometric  stiffness  matrix.  The  updating 
of  the  basic  matrix  and  the  geometric  matrix  are  not  necessarily  done  at  the 
same  time,  because  the  geometric  matrix  contains  the  stresses  themselves.  It 
is  not  satisfactory  to  update  the  geometric  stiffness  matrix  when  the  stresses 
have  relatively  large  errors.  For  this  reason,  the  updating  of  the  geometric 
stiffness  matrix  is  only  permitted  when  the  error  state  of  the  solution  is 
within  certain  bounds  controlled  by  parameters  within  the  code. 

It  was  felt  worthwhile  to  evaluate  the  new  nonlinear  element  in  comparison 
with  conventional  finite  element  approaches,  In  application  to  nonlinear 
analysis.  In  order  for  such  a  comparison  to  be  valid,  it  must  be  made  with  the 
two  element  types  having  all  features  In  common,  i.e.,  solution  procedure, 
element  nodal  and  freedom  formulation,  all  updating  and  transformations, 
etc.,  except  for  the  single  basic  feature  which  distinguishes  the  new  element: 
the  use  of  higher  order  axial/membrane  freedoms  in  combination  with  lower 
order  bending  freedoms.  In  order  to  accomplish  the  required  comparison,  the 
new  element  was  provided  with  a  special  solution  procedure  option:  a  solution 
procedure  constraint  was  developed  which  constrains  the  axial  freedoms  at 
nodes  2,  3  and  4  to  values  which  are  the  Interpolated  values  at  these  nodes  of 


14 


? 


% 


an  element  *hose  axial  displacements  are  linear  functions  (a  conventional 
axial  displacement  function).  This  constraint  process  is  of  the  type  called 
"multi -point  constraint",  and  Is  implemented  by  a  set  of  matrix  transforma¬ 
tions.  It  permits  the  element  to  be  used  in  either  its  "nonlinear"  mode  of 
behavior,  or  in  a  conventional  mode,  with  all  other  computational  processes 
unchanged.  As  noted  in  Table  1,  this  option  is  effected  by  input  of  the  value 
0PT2  *  negative.  Calculations  made  with  this  option  converged  very  slowly  to 
erroneous  results,  verifying  that  the  conventional  type  of  element  is 
unsuited  to  the  accurate  computation  of  highly  nonlinear  problems.  Section 
2.3  discusses  these  results,  which  are  in  complete  agreement  with  the  earlier 
results  of  Haftka,  Mallett  and  Nachbar  (Reference  1).  This  capability  was 
used  in  the  static  perturbation  version  of  the  code,  and  is  not  currently 
operative  in  the  BFGS  version. 

Theory  was  developed  for  a  beam  element  capable  of  bending  and  twisting  In 
three-dimensions,  called  herein  the  3-0  beam  element.  This  type  of  element  is 
required  for  problems  such  as  the  nonlinear  bending  and  torsional  deformation 
and  buckling  of  beam  structures.  A  difficult  technical  problem  arises  at  the 
outset  of  this  type  of  derivation.  It  is  recalled  that  one  of  the  goals  of  the 
present  research  is  to  avoid  unrecoverable  cunulative  error  in  the  problem 
solution.  A  principle  offender  in  this  regard  is  the  calculation  of  the 
strain  itself.  To  avoid  such  errors,  the  deformation  must  be  determined  by 
direct  calculation  of  the  total  strain,  using  the  total  displaced  state  of  the 
structure,  rather  than  by  strain  incrementation.  To  do  this  requires  a 
precise  definition  of  the  rotation  state.  The  three-dimensional  beam  element 
undergoes  three  components  of  rotation.  These  include  the  twist  and  the  two 
bending  rotations,  all  of  which  can  be  large  for  nonlinear  problems.  In  the 
large  rotation  state  the  orientation  of  a  rotated  element  or  a  node  of  an 
element  in  space  cannot  be  represented  by  arbitrarily  ordered  cartesian  com¬ 
ponents  referred  to  a  fixed  coordinate  system.  Neither  can  the  orientation  be 
arrived  at  by  sumning  small  rotations  referred  to  cartesian  systems.  The 
basic  problem  is  that  rotations  are  not  vectors  and  therefore  are  not  additive. 


15 


A  correct  large  rotation  state  can  be  obtained  through  the  use  of  sequenced 
angular  rotations  called  Euler  angles.  Each  Euler  rotation  takes  place  about 
an  axis  which  has  been  subjected  to  all  prior  rotations  In  the  sequence,  which 
must  take  place  In  a  specified  order.  This  approach  has  been  successfully 
used  for  the  calculation  of  large  motions  of  spacecraft  as  well  as  other  types 
of  large  rotation  dynamic  problems  (Reference  9).  It  Is  necessary  to  use  this 
approach  to  develop  the  3-0  beam  element.  If  the  conventional  small  angle 
(cartesian)  approach  were  used,  the  strains  would  be  inconsistent  with  the 
rotation  state  of  the  structure,  and  it  would  be  Impossible  to  correct  the 
equilibrium  configuration  through  the  residual  load  iteration  process. 

Strain  displacement  relations  based  on  Euler  angles  were  not  found  In  the 
literature,  and  consequently  a  set  of  appropriate  finite  element  deformation 
equations  had  to  be  developed.  The  derivation  was  carried  out  using  a  tensor- 
lal  approach  and  convected  coordinate  systems.  Appendix  B  describes  the 
approach  In  some  detail.  This  development  presented  a  number  of  difficulties, 
including:  the  need  develop  rational  approximations  associated  with  the 
relative  importance,  of  ..  many  different  types  of  nonlinear  terms  in  the 
strain  displacement  relations;  the  deformation-following  beam  cross-sectional 
axis  system  (Euler- angle- defined)  does  not  maintain  its  "longitudinal"  axis 
along  the  beam  centerline  axis,  and  It  is  necessary  to  define  an  additional 
Euler- angle- defined  convected  system  which  has  this  desirable  property  (see 
Appendix  B,  Figure  B5);  It  Is  necessary  to  determine  incremental  cartesian 
bending  and  twisting  angles  as  well  as  Euler  increments,  in  order  to  maintain 
physical  reality  in  Interpretation  of  the  deformations  and  the  loads.  In 
Implementing  the  3-D  element  together  with  a  stepping  solution  procedure  it  is 
necessary  to  use  a  large  number  of  transformations  of  geometrical  types,  in 
order  to  maintain  and  update  tne  different  convected  coordinate  systems  and 
the  two  sets  of  Euler  angle  totals.  One  such  transformation  provides  the 
needed  relationship  between  the  stepwise  Euler  angle  rotation  Increments  and 
the  cartesian  increments;  the  required  transformation  Is  called  the  if  trans¬ 
formation  (Reference  9).  Other  transformations  are  required  to  transform  the 
stiffness  matrix  of  the  element  from  Its  derivation  coordinate  system  to  the 
coordinate  systems  of  the  solution  process;  that  is,  to  transform  the  stiff¬ 
ness  matrix  from  those  definitions  used  in  the  strain  displacement  relations 
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to  those  which  are  suitable  for  merging  together  adjacent  elements  and  obtain¬ 
ing  the  problem  solution  In  terms  of  meaningful  cartesian  quantities.  Also, 

In  the  solution  procedure  the  coordinate  systems  used  for  each  element  are 
convected,  that  Is,  they  follow  the  elements  throughout  the  deformation  proc¬ 
ess  In  order  to  retain  for  each  elanent  a  small  deformation  state.  Hence,  It 
Is  necessary  to  perform  repeated  transformations  to  accomplish  the  geometri¬ 
cal  updating  of  various  solution  and  geometry  parameters.  The  total  solution 
process  for  the  3-0  beam  element  has  been  flow-charted  to  provide  a  methodol¬ 
ogy  description  suitable  for  computer  coding.  A  computer  code  for  this 
element  has  been  about  75%  completed. 

2.3  Illustrative  Nimerlcal  Results 

This  section  presents  nimerlcal  results  for  several  beam  bending  problens. 
Results  obtained  with  both  the  static  perturbation  method  and  the  BFGS  method 
are  discussed.  The  purpose  of  the  example  problems  is  to  Illustrate  the 
displacement  magnitude  capabilities,  convergence  characteristics,  and  limita¬ 
tions  of  the  methods  developed.  Refer  first  to  Figure  1.  A  simple  two 
element  beam  structure  Is  bent  by  an  end  load  In  the  global  Z  direction.  The 
loaded  end  of  the  beam  is  either  completely  free  (Figure  la)  or  constrained 
against  X  -direction  displacement  (Figure  lb).  This  problem  was  solved  by  the 
static  perturbation  approach. 

Tables  2-6  present  numerical  data  for  the  beam  structure  of  Figures  la  and  lb. 
The  tables  Include  deflections,  rotations,  axial  forces  in  the  elements,  and 
convergence  data  (numbers  of  Iterations  and  percent  error  based  on  residual 
loads).  The  two  values  of  axial  force  shown  In  the  tablas  for  node  5  are  those 
computed  for  the  two  elements  which  connect  to  this  mode. 

Table  2  gives  nimerlcal  results  for  the  problem  of  Figure  la.  The  load  varies 
from  0  to  720  (pounds),  while  the  end  deflection  varies  from  0  to  2.70 
(Inches).  A  graph  of  the  deflection  versus  the  load  would  be  very  nearly  a 
straight  line,  as  the  only  nonlinearity  In  this  problem  results  from  the  small 
foreshortening  of  the  beam  due  to  its  deflection  and  rotation.  The  deflection 
of  2.70  In.  Is  93%  of  the  theoretical  value  for  this  beam,  a  reasonable  value 
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for  a  two  element  model  where  each  element  has  only  a  second  degree  displace¬ 
ment  function  capability  In  bending.  The  average  rotation  of  element  #2  at 
720  lbs.  is  .189  radians.  The  average  element  rotation  is  the  angle  to  which 
the  "baseplane"  of  the  element  Is  updated  In  geometrical  updating.  The  table 
shows  the  sequence  of  baseplane  update  angles  for  both  elements.  The  deflec¬ 
tion  of  2.70  In.  Is  13.5%  of  the  total  beam  length,  an  amplitude  which  Is  well 
Into  the  potentially  nonlinear  range  in  structural  analysis.  However,  since 
the  right  end  of  the  beam  Is  permitted  to  deflect  freely,  the  beam  is  not 
stressed  axially  due  to  nonlinear  strain  buildup,  and  the  behavior  is  essen¬ 
tially  linear.  The  payoff  of  the  new  nonlinear  element  in  this  case  Is  that 
It  al  1  ows  the  axial  stress  to  ignore  the  nonlinear  strain  effects,  even  though 
fully  nonlinear  strain  calculation  is  done  in  the  analysis.  This  is  accom¬ 
plished  by  the  quartic  ax*‘al  displacement  function. 

The  element  axial  stresses  are  small  and  essentially  constant  over  each  ele¬ 
ment.  The  jump  in  axial  stress  at  node  5  balances  a  corresponding  jump  in  the 
shear  at  this  node.  Nearly  constant  axial  stress  In  an  element  is  required  by 
the  equilibrium  equations.  The  nonzero  axial  stress  values  are  correct,  and 
result  from  the  inclination  of  the  end  of  the  beam  with  respect  to  the  applied 
load.  This  is  Illustrated  in  Figure  2.  The  figure  gives  the  equations  of 
equlllbrlur  h  ilch  mu-t  be  satisfied  by  the  shear  and  axial  forces  at  the  end 
of  the  beam.  It  is  seen  that  the  inclination  of  the  shear  force  requires  the 
axial  stress  4n  the  beam  to  be  nonzero.  The  inclination  of  the  shear  force 
can  be  accounted  for  either  in  the  element  strain  formation  or  In  the  solution 
procedure  (by  geometrical  uodating).  To  account  for  this  in  the  strain 
formulatior,  it.  Is  necessary  to  retain  nonlinear  tens  in  the  beam  transverse 
shear  strain,  a  type  of  nonlinearity  not  usually  retained  In  geometrically 
nonlinei,**  analysis.  A  simpler  approach,  and  that  used  to  compute  the  data 
under  discussion  here,  is  to  use  the  geometrical  updating  to  rotate  the  shear 
force.  The  result  Is  that,  at  the  72(1  lb  load,  707  lbs.  is  normal  to  the  beam 
(the  shear,  S)  and  135  lbs  Is  directed  along  the  beam  axis,  producing  the 
axial  stress  shown  in  Table  2.  These  values  are  determined  by  the  inclination 
of  the  updated  baseplane  which  Is  seen  In  the  table  to  be  .189  radians.  The 
total  load  in  the  global  Z  direction  remains  720  lbs.  It  is  noted  that  at 
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larger  displacements  the  element  stresses  do  not  necessari ly  fol low  the  sim¬ 
ple  line  of  explanation  given.  For  such  problems  the  solution  converges  to 
give  apparently  erroneous  stresses  in  some  cases. 

Faster  convergence  can  usually  be  obtained  if  geometrical  updating  is  done 
Infrequently.  This  Is  because  the  curved,  updated  geometry  of  the  beam  causes 
the  axial  and  bev'*ng  displacements  to  be  numerically  "coupled"  much  more  than 
they  are  in  the  initial  flat,  nonupdated  geometry,  and  therefore  slows  the 
Iterative  con^rgence  process.  TMs  suggests  that  solutions  might  be 
obtained  at  lower  cost  by  applying  the  total  load  In  the  first  load  step.  In 
this  approach  only  one  geometrical  updating  Is  required,  corresponding  to  the 
final  deflected  state.  Such  a  solution  is  shown  on  Table  3.  The  element  #2 
baseplane  was  updated  in  one  step  to  the  Inclination  of  .193  radians,  con¬ 
sistent  with  the  zero^-  Iteration  displacements,  for  which  the  displacement 
at  the  end  of  the  beam  is  2.75  in.  (final  convergence  was  to  2.69  in.).  The 
computed  displacements  are  almost  identical  to  those  of  Table  2.  The  axial 
stresses  are  slightly  different  because  the  baseplane  Is  updated  to  a  slightly 
different  angle  than  the  .189  radians  of  Table  2.  The  results  of  Table  3  show 
two  Important  facts:  the  converged  result  for  large  loads  can  be  obtained  In 
a  single  load  step;  certain  aspects  of  the  solution,  such  as  the  axial  stress 
in  the  present  problem,  may  be  sensitive  to  the  inclination  of  the  updated 
baseplane,  so  that  updating  should  not  In  general  be  neglected. 

Geometrical  updating  Is  o.ily  required  when  the  element  rotation  relative  to 
the  baseplane  coordinate  systems  become  large,  e.g.,  greater  than  about  20°, 
Table  4  shows  a  case  of  delayed  updating.  Here  the  updating  has  been  delayed 
until  the  load  reaches  660  lbs.  The  deflection  results  are  nearly  the  same  as 
those  of  Table  2,  but  of  course  the  axial  stresses  are  not  correctly  computed 
until  the  baseplane  Is  updated.  Note  that  the  final  baseplane  angles  here 
have  resulted  from  a  less  recent  update,  and  hence  differ  slightly  from  those 
o'  Table  3. 

The  problem  of  Figure  lb  has  the  global  X  displacement  constrained  at  the 
right  end.  This  problem  Is  highly  nonlinear,  typical  In  character  to  many 
practically  Important  cases  Involving  end-or-edge-constrai ned  beams  and 
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plates.  The  solution  data  are  given  In  Table  5  and  on  Figure  3  for  a  load 
range  of  0  to  240  lbs.  The  baseplane  In  this  example  Is  updated  when  the 
average  slope  of  the  beam  (relative  to  the  most  recent  baseplane  update) 
exceeds  .01  radians.  Convergence  Is  reasonably  fast  except  at  the  180  lb. 
load.  The  axial  stresses  are  essentially  constant  over  the  lengths  of  the 
elements,  being  dominated  by  the  effect  of  the  end  constraint.  They  are 
responsible  for  the  nonlinear  stiffening  behavior  illustrated  by  the  force- 
deflection  plot  of  Figure  3.  For  larger  load  levels,  the  degree  of  non¬ 
linearity  of  this  problem  Increases  very  rapidly. 

The  next  example  concerns  a  "conventional"  beam  element.  This  element  Is 

Identical  to  the  new  nonlinear  element  except  that  the  axial  displacements  at 

the  Interior  nodes  are  constrained  to  take  values  defined  by  linear  Interpo- 
« 

latlon  between  the  end  node  values.  That  is,  these  freedoms  are  In  effect 
omitted  from  the  problem  solution  by  constraining  the  axial  displacement 
shape  to  be  the  linear  shape  of  the  "simplex"  type  of  element.  In  solving 
problems  with  the  constrained  element,  fully  nonlinear  axial  strains  due  to 
the  bending  displacements  are  retained,  and  the  new  nonlinear  solution  pro¬ 
cedure  Is  also  used.  Thus,  the  results  provide  a  consistent  comparison 
between  the  new  type  of  element  and  one  of  conventional  formulation,  with  all 
other  aspects  of  the  numerical  processing  kept  the  same.  The  results  are 
tabulated  In  Table  6  for  the  load  range  of  0  to  240  lbs.  and  plotted  on  Figure 
4.  The  solution  at  720  lbs.  was  computed  In  a  single  load  step.  The  final 
deflection  at  720  lbs.  Is  1.57  In.,  which  Is  considerably  less  than  the  2.70 
in.  of  the  new  nonlinear  element.  This  error  reflects  the  excessive  and 
erroneous  stiffness  of  the  conventional  element  due  to  the  axial  stresses 
which  are  "locked"  In  this  element  by  nonlinearity.  The  "simplex"  axial 
displacements  cannot  ronove  these  locked-in  axial  stresses  because  of  defi¬ 
ciencies  of  their  functional  forms.  The  axial  stresses  are  seen  in  the  table 
to  be  very  large.  The  error  Illustrated  by  this  example  is  consistent  with 
that  discussed  in  Reference  1.  The  new  nonlinear  element  has  eliminated  this 
type  of  error. 

For  displacements  and  rotations  which  are  significantly  larger  than  those  of 
the  previous  examples,  convergence  difficulties  were  encountered  with  the 
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static  perturbation  approach.  As  discussed  In  Section  2.2,  the  BFGS  method 
was  implemented  to  Improve  convergence.  In  addition,  the  transverse  shear 
strain  was  modified  to  Incorporate  nonlinear  terms.  The  primary  effect  of 
this  modification  Is  to  rotate  the  resultant  shrar  force  on  the  element  so 
that  It  is  parallel  to  the  deformed  beam  cross-section.  This  nas  a  relatively 
small  effect  on  problem  solutions. 

Figure  5  shows  two  cantilever  beam  problems  which  were  solved  with  the  PCGS 
approach.  Figures  5a  and  5b  show  a  single  element  problem,  with  the  support 
located  both  at  the  center  of  the  element  and  at  the  left  end.  The  problem  of 
Figure  5a  requires  no  geometrical  updating  because  the  baseplane  does  not 
rotate.  It  yields  an  exact  solution  for  the  rotations  of  the  ends  of  the 
beam.  In  contrast,  the  problem  of  Figure  5b  has  significant  element  baseplane 
rotations  requiring  geometrical  updating.  The  loading  in  both  cases  is  a  pure 

moment,  and  the  purposes  of  the  example  are  to  investigate  the  rate  and  degree 

of  convergence  obtainable  at  large  displacement  and  the  influence  of  geo¬ 
metrical  updating  on  convergence  and  accuracy.  Figure  6  shows  dimensioned 
sketches  of  the  deflections  for  both  cases.  For  each  problem,  the  figure  also 
shows  the  displaced  condition  referred  to  the  coordinate  system  of  the  other 
problem.  These  data  are  shown  in  parentheses.  The  deflections  are  large,  on 
the  order  of  half  the  length  of  the  beam.  It  is  seen  that  the  simple 

cantilever  element  converges  to  a  solution  having  4%  to  5%  more  curvature  than 

that  for  the  doubly  cantilevered  (synmetrical )  element.  The  cause  of  the 
difference  is  almost  certainly  the  geometrical  updating  required  for  the 
simple  cantilever  case,  which  Is  not  done  in  the  symmetrical  case.  Conver¬ 
gence  for  the  symmetrical  case  occurs  on  the  first  iteration;  it  is  much 
slower  for  the  case  with  geometrical  updating.  The  values  of  the  element 
shear  and  axial  stresses  for  the  simple  cantilever  differ  from  those  of  the 
symmetrical  problem.  This  causes  element  (residual)  loads  which  Influence 
solution  convergence  and  accuracy,  particularly  for  multi-element  problems. 

The  cause  of  the  numerical  differences  between  these  two  solutions  has  not 
been  fully  resolved.  It  appears  doubtful  that  the  differences  arise  due  to 
the  element  formulation  Itself,  because  the  symmetrical  problem  yields  an 
exact  solution.  The  simple  cantilever  case  can  clearly  have  exactly  the  same 


21 


deformation  state  as  the  double  cantilever,  when  referred  to  an  updated  base- 
plane.  Thus,  the  determination  of  residual  loads,  which  is  based  on  only 
displacements  referred  to  the  updated  baseplane,  is  potentially  Identical  for 
the  two  problems.  It  is  likely  that  the  BFGS  solution  procedure  (which  Is 
basically  an  optimum-seeking  type)  has  become  trapped  along  a  solution  path 
which  has  a  false  minimum  error  state.  This  view  Is  strengthened  by  the  fact 
that  frequently  a  regeneration  cf  the  stiffness  matrix  in  the  deformed  state, 
followed  by  subsequent  BFGS  calculations,  yields  a  substantially  improved 
solution  accuracy.  The  extensive  use  of  single  precision  arithmetic  in  the 
code  may  also  be  a  factor  in  the  minimum  error  problem.  It  Is  also  noted  that 
the  geometrical  updating  of  the  axial  displacement  values  at  nodes  2  and  4 
makes  use  of  interposed  bending  displacements  at  these  nodes.  The  bending 
displacements  themselves  are  not  updated  for  these  nodes.  However,  if  they 
were  updated,  the  resulting  values  would  not  adhere  to  the  quadratic  bending 
displacement  form  when  referred  to  the  updated  baseplane.  This  is  a  source  of 
Inconsistency  inherent  In  the  uodating  process,  due  to  the  different  function 
shapes  used  for  the  axial  and  bending  displacements.  This  inconsistency 
should,  however,  be  correctable  by  residual  load  iteration. 

2.4  Suggested  Further  Work 

The  element  development  appears  to  have  successfully  controlled  the  problem 
of  large  nonlinear  strains,  as  illustrated  by  the  results  shown  on  Figure  6 
for  the  symmetrical  cantilever.  However,  for  more  general  cases,  such  as  that 
of  the  unsymmetrlcal  case  on  Figure  6  and  a  number  of  multi-element  problems 
which  have  been  solved,  there  remain  unresolved  problems  in  either  the  formu¬ 
lations  or  In  the  actual  calculations.  Displacements  which  are  very  large 
have  been  successfully  computed  despite  the  difficulties  encountered,  how¬ 
ever,  and  the  potential  of  the  new  type  of  element  appears  to  have  been 
adequately  demonstrated.  If  the  research  is  to  be  continued,  the  accuracy  and 
convergence  difficulties  will  have  to  be  addressed  as  a  first  step.  The  four 
areas  of  study  listed  below  are  suggested  candidates  for  this  work. 


o  eliminate  the  use  of  single  precision  arithmetic  in  all  suspect 
calculations. 


o  investigate  the  possibility  of  the  solution  procedure  being 
"trapped"  by  a  false  minimum  of  the  error  level. 

o  investigate  whether  the  limitation  of  the  element  to  constant  cur¬ 
vature  (w  is  quadratic)  while  simultaneously  a  quadratic  rotation 
is  allowed  is  a  cause  of  inconsistency  which  could  contribute  to 
numerical  problems. 

It  is  expected  that  the  numerical  difficulties  which  still  exist  can  be 
resolved.  In  this  event  it  appears  particularly  important  to  develop  and  test 
the  three-dimensional  beam  element.  This  work  will  evaluate  the  Euler  angle 
deformation  theory  and  the  set  of  geometrical  transformations  inherent  in 
this  approach  which  is  new  in  the  field  of  nonlinear  finite  element  stress 
analysis.  The  theory  and  procedural  specifications  have  been  completed  and 
the  computer  coding  partially  completed  for  this  task. 

If  the  three-dimensional  beam  element  work  shows  the  Euler  angle  theory  to  be 
a  valuable  tool  for  large  displacement  finite  element  analysis,  consideration 
should  be  given  to  a  further  task.  This  task  would  apply  the  Euler  angle 
approach  and  the  algorithms  developed  for  the  two-dimensional  beam  element  to 
the  development  of  plate  and  shell  elements.  It  would  result  in  a  truly  large 
deflection  analysis  method  for  plates  and  shells,  leading  eventually  to  a  much 
needed  large  displacement  shell  buckling  analysis  capability. 
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Figure  1:  Illustrative  Cantilever  Beam  Problems 
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Figure  2:  Beam  End  Conditions  for  Large  Rotations 


LOAD,  Pounds 

Deflection  Behavior  for  Restrained  -  End  Beam  (Figure  lb) 
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Figure  4:  Comparison  of  Load-Deflection  Behavior  for  Nonlinear  and  Conventional  Elements 
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Figure  6a:  Simple  Cantilever 
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for  Problem  of  Figure  6a 
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Figure  6b:  Double  (Symmetrical)  Cantilever 


Figure  6:  Displacements  of  Single  Element  Cantilever  Beam 
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Table  1:  Summary  of  Nonlinear  Solution  Procedure  Controls 
(a)  OPTION  CONTROLS  0> 


STATIC 

PERTURBATION 

ORDER 

PATH  PARAMETER 

"ALTERNATE" 

ITERATION 

LINEARIZED 

SOLUTION^ 

"CONVENTIONAL" 

ELEMENT 

LOAD 

TYPE 

~ 

EH 

Kail 

NO 

YES 

1  B> 

0PT2-1 

0PT1--2 

No 

OPTI* 

-1 

OPTl- 

1 

Yes 

0PT2*0 

Yes 

0PT2--1 

2 

Yes 

0PT2-2O 

Yes 

0PT2-2 

0PT1* 

-1 

OPTI* 

1 

No 

No 

3 

Yes 

0PT2-3O 

Yes 

0PT2-3 

OPTI* 

-1 

OPTI* 

1 

No 

No 

E>  OPT!  and  0PT2  are  program  Input  data. 

It>  First  order  static  perturbation  Is  standard  stepwise  solution  - 
Path  parameter  Is  always  load-type  In  this  case 

B>  Linearized  solution  omits  all  types  of  nonlinearities. 

E>  Partial  list  given;  Includes  principal  control  inputs. 


(b)  NUMERICAL  CONTROLS 


INPUT  DATA 
NAME 


PROCESS  CONTROLLED 


BY  INPUT  DATA 


EPSCON; 

Allowable  Error  (Residual  Loads)  For  Convergence 

DALIM: 

Rotation  of  Element  From  Base  Plane  at  Which 

Geometrical  Update  is  Performed. 

EQCHK: 

Rate  of  Divergence  at  Which  Stiffness  Matrix  Update 
is  Performed. 

Error  (Residual  Loads)  at  Which  Stiffness  Matrix  Update 
is  Performed. 

ITRLIM,,  Iteration  and  Updating  Counts  at  Which 
UPLIM  •  Solution  is  Aborted _ _ 


Required  Minimum  Number  of  Iterations  Between  StiffnesT 
Matrix  Updates  (Over-Rides  EQCHK.  RACHK) _ 
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APPENDIX  A 


Static  Perturbation  Path  Parameter 


This  Appendix  outlines  briefly  the  use  of  the  load-based  and  deformation- 
based  path  parameters  in  the  static  perturbation  method.  The  equations  given 
herein  are  programmed  in  the  2-D  beam  code  (see  Section  2.2). 

The  static  perturbation  method  uses  Taylor  series  to  represent  the  incre¬ 
mental  displacement  vector  aq  and  the  incremental  load  vector  aP  . 

fcP=  PS  +  +  jp  s  + 

.  ...  i.  ,  •“  3 

CbQ  *  QS  +  -jQS  +  -  Q  S  +  •  *  • 


where  (  *  )  denotes  differentiation  with  respect  to  the  path  parameter  S.  It 
is  convenient  to  represent  the  P  derivatives  as  follows: 

•  •  A 

P  «  -X  P 

•«  ••  e 

psi.? 

•  ••  a 

P  *  A.  P 


so  that 


A.P  =  p’(X*  +  1-^-^  +  £  X  s* 
and  to  set  P°  equal  to  the  load  increment 
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such  that 


«*  i  ^ 

a  s  ♦  4  -A  *  =  t 

The  values  of  the  Q  derivatives  can  be  shown  to  be  given  by 

•  -<  !  a* 

<5  a  K  J\.f 

a=  k'( Ip*  -  *<*) 


•  •• 

in  which  K  is  the  tangent  stiffness  matrix  and  K  and  K  are,  respectively,  the 
first  and  second  derivatives  of  the  K  matrix  with  respect  to  the  path  param¬ 
eters  (accomplished  by  chain-rule  differentiation:  >k  AS  q) . 

The  solution  process  requires  solution  first  for  Q,  followed  by  calculation  of 
‘K,  followed  by  calculation  of,  in  order,  Q,  K,  and  finally,  Q.  The  Taylor 
series  then  gives  the  value  of  the  vector  &.Q. 

The  load-based  path  parameter  sets 

JS_  «  A.  *  O  f  \ 

6,p«  P*JuS  s  5  P* 

S-l 

«• 

and  allows  calculations  of  &Q  immediately  that  Q  and  Q  are  known.  The  first, 

2 

second,  and  third  order  approaches  retain,  respectively,  the  terms  S,  S  ,  and 
$3.  This  is  a  relatively  simple  approach  to  implement. 

The  deformation  -  based  path  parameter  is  much  more  complex  because  the  values 
of  JU.A,  and  S  are  unknown.  This  approach  is  based  on  the  definition  of  S 
below,  using  the  tangent  stiffness  matrix, 

e  fcCl-  K  *0. 
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where  (  )T  denotes  the  row-vector  (transpose).  Thus,  S  Is  roughly  propor¬ 
tional  to  the  incremental  displacement  amplitude.  To  explain  the  solution 
process  for  this  case  it  Is  necessary  to  make  a  number  of  definitions,  as 
given  below. 

O  I 

Q  »  V<-  ? 

(evaluate  K  with  Q°  In  place  of  Q) 

o>*=  k°  ct  ;  x'pl- 

(evaluate  K  with  Q°  in  place  of  Q  in  the  Q  - 
dependent  part  of  K,  called  Kj) 

(evaluate  K  with  Q°  in  place  of  Q  in  the  Q  - 
dependent  part  of  K,  called  here  K^) 

(evaluate  K  with  Ql°  in  place  of  Q  in  K^l 


Rf 

kV- 

K3°  - 


with  these  definitions, 

•»  ,1 

K  *=  -A_ 


*’  *a  kV  +XkV  ••  <V 


and  (see  equation  for  Q),  x  t 

+  k  tk<a  «  a. A.K*  (-A.Q.  -  -A-  Q.'  )  +  x 

JLa*  (i?  k\*  -h  X  kV  -  X.  k*°  ) 
It  is  convenient  to  rewrite  ohis  as 

,K,a*E»i-  -0.X. (»"*-* '**)  *  -0s (fi’-m'-fO 


where  the  various  "PN°"  are  derived  by  substitutions  for  the  various  "K"  type 
matrices  above. 

Finally,  we  define  the  vectors 

(ipi#  pV) 

«OQ*  a  |C‘(P3“-  2-PZ°  ~  P5*) 
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.»is  il  n  i  » -  i 


with  the  result 

k  q  -  Ik.  at  -  _a_  (^a.a°  —  .A--A-  q'<s- 

•  M  »•* 

It  is  seen  that  the  entire  evaluation  depends  on  knowi ng A  and  S. 

These  values  can  be  computed  (with  much  difficulty)  from  the  basic  definition 
of  S2  given  above.  The  equations  are  complicated  are  are  omitted  here.  The 

final  result  is  ,  % 

=  of  -  ai*  ~  aa"  U  JtJVS  j  -aaa  5  ) 


The  second  order  procedure  keeps  only  S2,  end  the  third  order  procedure  keeps 
both  S3  terms. 
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Euler  Angle  Theory  for  Beam  Element 

This  Appendix  describes  the  use  of  Euler  angles  in  the  determination  of  the 
deformations  of  beam  elements,  emphasizing  the  physical  nature  and  basis  for 
this  approach.  The  derivation  of  the  strain-displacement  equations  in  terms 
of  the  Euler  angles  is  also  described  briefly. 


Figure  B1  shows  a  beam  cross-section  in  the  initial  undeformed  state.  The 
section  is  shown  rectangular  only  to  aid  visual  clarity.  The  xyz  triad  is 
oriented  such  that  x  is  the  beam  centerline  and  the  y  and  z  axes  are  the  axes 
of  bending.  The  shear  deformations  associated  with  bending  occur  in  the  xy 
and  xz  planes. 


Figure  B2  shows  the  cross-section  after  displacement.  The  triad  xyz  has 
followed  the  motion  of  the  section,  as  described  below.  This  "convected"  xyz 
triad  is  orthogonal,  but,  as  will  be  shown  below,  is  not  truly  "imbedded"  in 
the  material.  The  xyz  triad  is  carried  Into  the  xyz  triad  by  means  of  the 
sequence  of  Euler  angle  rotations,  in  the  following  manner. 

Allow  first  a  rotation  /3 ^  of  the  undeformed  section  about  the  x  axis,  as 
shown  in  Figure  B3a.  This  results  in  a  new  triad,  denoted  in  the  figure  by 
x'y'z'.  Of  course,  x'  and  x  are  identical. 
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(a)  Rotation  /?,  About 
x  Axis 


(b)  Rotation  0,  About 
y'  Axis 

(c)  Rotation  /?,  About 
z"  Axis  J 


Figure  83  -  Euler  Rotations 


Next  a  rotation  0^  takes  place  about  the  y'  axis,  resulting  in  a  new  triad 
x"y"z".  Finally,  a  rotation  f&^  about  z"  results  in  che  triad  xyz  which 
describes  the  deformed  orientation  of  the  cross-section.  The  angles  (5^, 
and  are  Euler  angles.  They  are  restricted  to  the  sequence  of  axes  (1-2-3) 
in  the  special  meaning  and  sequence  illustrated.  For  the  general  case,  there 
are  12  possible  Euler  angle  sequences,  but  only  the  (1-2-3)  sequence  described 
above  is  needed  for  the  present  discussion.  The  importance  of  this  rotation 
description  in  the  present  application  is  that  it:  (1)  fully  accounts  for  the 
effect  of  large  rotations  in  reorienting  the  beam  cross-sectional  axes;  (2) 
avoids  any  errors  due  to  cartesian  rotation  incrementation. 

The  three  Euler  rotations  are  not  in  general  those  of  conventional  beam  twist 
and  bending,  although  for  small  total  rotation  magnitudes  they  are  indistin¬ 
guishable  from  these  quantities.  For  large  rotations  they  are  not  correctly 
viewed  in  this  way,  however,  and  for  this  reason  it  is  incorrect  to  attribute 
the  beam  twisting  and  bending  stif .  ness  properties  to  the  rotation  values  (3  y 
(1  2,  and  Derivation  of  the  beam  twist  and  bending  moments  instead  are 

derived  by  a  rigorous  process  described  later.  While  the  Euler  angle  repre¬ 
sentation  has  the  advantage  of  rigor,  it  does  not  provide  a  fully  satisfactory 
deformation  description  from  a  physical  viewpoint.  To  obtain  the  needed 
physical  interpretation,  it  is  necessary  to  compute  small  incremental  rota¬ 
tions  about  the  beam  bending  and  twisting  axes,  superimposed  on  a  previously 
accumulated  large  rotation  state.  Figure  B4  attempts  to  illustrate  this  view 
of  the  deformation.  The  figure  shows  the  deformed  section  with  the  associated 
triad  ixyz.  The  triad  is  essentially  identical  with  the  beam  section  axes  and 
centerline,  deviating  only  slightly  from  these  axes  due  to  the  shear  strains 
(angles  of  the  order  of  0.3  degrees).  The  figure  indicates  that  a  sma'll 
rotation  superimposed  on  the  section  in  its  xyz  orientation  can  be  viewed  in 
two  ways:  as  a  cartesian  rotation  taking  place  about  the  x,y,  and  z  axes,  or 
as  an  increment  in  the  Euler  angle  values  (3  y  (3y  and  /3  y  The  Sfp 
i  fiy  arK*  S/2  3  be  viewed  as  taking  place  about  the  axis  systems  which 
were  defined  duri ng  the  entire  deformation  process-the  xyz  and  x'y'z'  and 
x"y"z"  systems.  This  is  indicated  in  the  figure.  These  rotations  are  not 
physically  meaningful,  but  can  be  used  in  a  mathematically  rigorous  way  to 
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P^  andd^S-j 


Figure  B4  -  Small  Superimposed  Rotations 


[<5/Jj  »  o^2 »  <5/^] 

[ <50$,  60a,  60a] 


describe  the  nonlinear  deformation  state.  On  the  other  hand,  the  cartesian 

small  rotations  SQ$,  £©-,  and  $6$  are  physically  meaningful  because 

x  y  z 

they  indicate  increments  of  twist  and  bencing  deformation;  and  in  dynamic 
analysis  they  are  correctly  associated  with  the  beam  cross-section  rotational 
inertia  properties.  To  combine  the  rigorous  nonlinear  deformation  description 
with  the  physically  meaningful  one  a  special  type  of  geometrical  transformation 
is  available.  This  is  expressed  by 

=  Wl \sfx  i 

W 


The  'TT  matrix  is  a  function  of  the  total  accumulated  and  fi A 

suitable  nonlinear  analysis  approach  must  be  based  on  deformations  described 
rigorously  by  0  and  i  p  ip  and  Addition  contain 

transformations  to  provide  numerical  results  in  the  form  £©$,  *©z* 

The  'IT  matrix  is  the  means  for  accomplishing  this. 

A  refinement  of  the  above  approach  cafl  be  made  to  eliminate  the  small  approxi¬ 
mation  in  the  meaning  of  the  ?©  values  which  is  due  to  the  fact  that  the  xyz 
triad  is  not  truly  identical  to  the  conventional  beam  twist  and  bending  axes. 


B4 


Figure  B5  shows  a  beam  cross-section  with  a  large  shear  deformation.  The 
actual  "material"  cross-section  is  shown  with  the  heavy 


Figure  B5  -  Beam  Cross-Section  Local  Coordinate  Systems 

lines.  The  light  dashed  lines  show  the  section  with  the  shear  deformation 
removed;  the  displacement  vectors  show  the  small  cross-section  motions  neces¬ 
sary  to  remove  the  shear  angles.  The  section  bisecting  axes  are  shown  for 
both  cases  as  an  aid  to  pictorial  clarity.  The  triad  xyz  coincides  with  the 
"mater ial-y  and  z"  cross-section  axes;  it  is  shown  by  the  very  solid  heavy 
coordinate  axis  lines.  The  x  axis  does  not  coincide  with  the  beam  material 
centerline.  The  triad  xyz*  coi nci  des  with  the  section  material-y  and  z  axes 
with  the  shear  deformation  removed.  Since  this  section  is  normal  to  the  beam 
centerline,  the  ~  axis  is  col  inear  with  the  centerline.  The  xyz  triad  is 
shown  by  the  very  heavy  dashed  lines.  It  is  desirable  to  keep  track  of  the  xyz 
triad  in  the  solution  process,  because  small  cartesian  rotation  increments 
referred  to  this  triad  are  precisely  the  "nominal"  section  twist  and  bending 
rotations.  This  is  easily  done  as  follows:  first  presume  that  Euler  angles 
^M1  ♦  ^M2’  ^M3  de^ine  the  triad  xyz;  then  (1)  compute  the  (see  fig.  B4) 
rotation  increment  S©*,  re-;;  (2)  compute  the  Euler  increment 

V  £^3  using  the  ir  matrix  and  sijn  to  obtain  the  total  p  p y 

Py  (3)  from  the  strain-displ acement  equations,  compute  the  shear  strain 
increment;  (4)  subtract  this  increment  from  the  Sep  and  se  j  values,  call¬ 
ing  the  results  and  66^  (these  are  the  bending  rotation  increments 
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which  correspond  to  conventional  beam  theory);  (5)  using  the  trans¬ 
formation  (the  IT  -type  transf ormation  which  is  computed  using  the  Euler 
angles  ^M2,  /^M3 ) »  compute  the  incremental  values  s { 

*^M3:  sun  to  compute  the  total  values  ^M2,  This  procedure 

maintains  the  totals  of  two  sets  of  Euler  angles:  /s ^2,  p  3  defining  the 
xyz  triad  (needed  to  compute  the  deformations);  and  /^M2*  P M3  definin9 

the  xyz  axes,  which  are  the  conventional  (convected)  axes  of  beam  theory. 

This  Appendix  closes  with  a  brief  outline  of  the  procedure  for  the  derivation 
of  the  nonlinear  strain-displacement  equations.  These  are  derived  in  terms  of 
the  Euler  angles  ^2,  /*3»  because  only  these  angles  provide  a  rigorous 
representation  of  the  total  rotations  of  the  beam  material  elements.  The  beam 
centerline  has  the  displacement  vector  <U-  .  Denoting  by  the  basis  vectors 
of  the  undeformed  beam  element,  and  by  x1  the  initial  undeformed  coordinates 
of  the  element,  the  initial  undeformed  beam  centerline  has  the  position  vec¬ 
tor,  respectively, 

-+■  VX  )  ©. ; 

The  vector  V  is  defined  to  be  the  additional  displacement  of  a  material 
point  which  is  off  the  centerline.  Hence  for  an  arbitrary  point  of  material 
on  a  beam  cross-section, 

=  (x:4  a')fti  -v  v 

The  Euler  angles  are  used  to  define  V  .  Referring  to  Figure  B3  and  defining 
the  basis  vector  sets  0^,  *nd  as  belonging,  respectively,  to 

triads  xyz,  x'y'z',  and  x"y"z",  v  can  be  written 


It  is  noted  that  actual  values  of  y'  and  y",  and  z'  and  z",  respectively,  for  a 
given  point  on  the  cross-section,  are  identical  to  the  values  of  y  and  z. 

This  is  because  all  of  these  coordinates  are  defined  by  "following"  the 
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material.  Hence,  with  y,  and  2  interpreted  as  initial  coordinates  of  a  point 
on  the  beam  cross-section,  V  is  more  simply  written 

V  »  1  *3 

The  values  of  and  «x  ^  can  be  represented  in  terms  of  the  initial  basis 
vectors  Cl^,  and  the  Euler  rotations  p  and  ^3.  The  derivation  is 

lengthy,  and  is  not  given  here.  The  result  is 

V*  /»,  (^3  -  1  °-0  + 

~  fj  c^(8»  +  «-i  »•*/*,  *••*>  ~  ®-3 

The  value  of  the  position  vector  R  is  now  known  for  any  point  on  the 
element,  in  terms  of  the  centerline  displacement  w.1 ,  the  Euler  rotations 
£3..,  and  the  "material"  coordinates  x,y,  and  z.  By  differentiating  R.  with 
respect  to  x,  y,  and  2,  there  are  obtained  an  important  set  of  vectors,  called 
the  basis  vectors  of  the  deformed  material  coordinate  system.  These  vectors 
contain  a  complete  description  of  the  deformation  state.  It  is  seen  that 
thesr  L-.-Js  vectors  contain  derivatives  of  the  p  i  and  the  ^  with  respect  to 
the  x,  y,  and  2  coordinates.  The  x-derivative,  in  particular,  is  important  in 
defining  the  deformation  of  the  beam  element.  The  theory  of  the  derivation 
process  follows  reference  13.  Simplifying  approximation  car.  made,  and  the 
details  and  results  are  too  lengthy  to  include  here.  It  is  simply  noted  that 
this  means  of  developing  the  strain  equations  is  exact  and  includes  all  of  the 
effects  of  nonlinearity. 

The  brief  description  of  the  nonlinear  beam  deformation  given  above  has  the 
purpose  of  illustrating  that  the  deformation  is  dete'mined  in  terms  of  the 
Euler  angles  rather  than  in  terms  of  conventional  cartesian  rotations.  This 
formulation  is  r  -ous  for  large  rotations  and  does  not  suffer  any  inaccu¬ 
racies  due  to  sunming  angular  motions. 
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APPENDIX  C 


Summary  of  Proposal  for  Contract  F49620-79-C-0057 

This  Appendix  contains  the  Technical  Approach  (Section  2)  and 
two  appendices  from  "Technical  Proposal  -  Program  for  Nonlinear 
Structural  Analysis",  submitted  to  AFOSR  In  August,  1978.  This 
proposal  Is  the  basis  of  the  current  AFOSR  contract  on  nonlinear 
finite  element  research.  The  discussions  herein  are  Intended 
to  provide  background  and  supplementary  Information  supporting 
our  present  report. 


2.0  TECHNICAL  APPROACH 


The  technical  approach  builds  on  existing  research  results.  The  element 
technology  to  be  used  is  basically  that  of  the  stability  elements  (Refer¬ 
ences  1,  2,  3).  The  solution  procedure  technology  will  be  based  on  the 
nonlinear-step  static  perturbation  procedure  (References  4,  5,  6).  The 
static  perturbation  procedure  has  been  demonstrated  to  be  a  superior 
solution  method  for  strongly  nonlinear  static  problems.  For  dynamic 
analysis,  an  extension  of  the  procedure  has  been  developed  in  the  current 
AFOSR  contract.  Numerical  data  demonstrating  the  superiority  of  this 
approach  are  given  in  this  proposal.  The  goal  of  this  proposed  research 
is  to  merge  these  two  technologies  into  a  working  pilot  computer  program. 

2.1  Technical  Requirements 

The  principal  features  required  of  the  overall  approach,  in  regard  to 
applicability  to  nonlinear  problems,  are  listed  below: 

Element  Technology  (References  1,  2,  3): 

1  Elements  are  required  whose  displacement  function  formula¬ 
tion  prevents  anomalous  (overstiff)  behavior  due  to 
nonlinear  strains.  These  are  called  stability  elements, 
and  utilize  extended  forms  of  axial /membrane  displacement 
functions,  in  conjunction  with  conventional  bending 
deformation  forms. 

2.  Element  strain  calculations  must  be  made  on  a  total 

strain  basis,  to  avoid  cumulative  errors  due  to  summing 
Increments.  This  is  required  to  allow  the  development  of 
large  rotations  and  nonlinearities. 
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3.  Element  displacement  functions  must  be  referred  to 
convected  coordinate  systems.  This  avoids  exchange  of 
axial /membrane  and  bending  displacement  roles  (e.g.,  u 
and  w  exchange  meaning  as  the  rotation  becomes  large)  in 
the  large  rotation  state,  and  permits  the  simplifying 
assumption  of  shallowness  in  forming  the  nonlinear  strain 
equations  for  shells  or  highly  deformed  beams  and  plates. 

4.  Residual  force  evaluation  and  equilibrium  corrections 
must  include  the  effects  of  element  strains  and  geometry 
changes. 

Static  Solution  Procedure  Technology  (References  4,  5,  6) 

1.  The  characteristic  problem  of  excessive  residual  forces, 
with  consequent  slow  convergence  or  divergence  in  problem 
solutions,  must  be  avoided  while  retaining  reasonably 
large  step  size.  This  requires  the  use  of  a  nonlinear 
stepwise  solution  procedure. 

2.  The  solution  procedure  must  be  compatible  with  the 
stability  elements,  in  particular  with  the  convected 
coordinate  system  approach. 

3.  The  solution  procedure  should  include  a  means  of  auto¬ 
matic,  internal,  computation  of  step  size.  Gains  in 
solution  economy  from  this  feature  can  be  very  large. 

Dynamic  Solution  Procedure  Technology 

1.  As  noted  above,  the  method  must  incorporate  a  nonlinear 
step.  Automatic,  internal  step  size  selection  should  be 
incorporated  insofar  as  is  possible. 

2.  The  solution  procedure  must  be  compatible  with  the 
stability  elements,  particularly  as  regards  the  convected 
coordinate  system  and  the  residual  load  calculations. 


C3 


•  • 


3.  The  solution  procedure  must  not  require  stepwise  Inver¬ 
sion  of  the  structural  stiffness  matrix.  Instead,  inver¬ 
sion  of  the  mass  matrix  must  be  used,  for  reasons  of 
economy. 

The  proposed  technical  approach  meets  all  of  these  goals.  All  of  the 
technical  developments  required  in  the  proposed  research  are  reasonably 
well  proven  as  regards  accuracy  and  practicability.  Hence,  their  merging 
Into  a  single  computer  program  appears  to  Involve  little  risk.  The 
major  gains  from  the  proposed  work  should  be  in  the  matter  of  evaluation 
of  the  overall  technical  approach  on  specialized  problems.  The  particu¬ 
lar  types  of  problems  for  which  this  approach  is  required  have  the 
following  characteristics: 

o  The  equilibrium  is  governed  primarily  by  nonlinear  axial /mem¬ 
brane  stresses  induced  by  bending  rotations. 

o  The  axial /membrane  stresses  vary  rapidly  over  the  structure.  . 
An  example  is  the  type  of  buckle  pattern  which  occurs  typi¬ 
cally  in  axially  compressed  cylinders,  in  higher  vibration 
modes  of  beams  and  plates,  and  in  short  wave  length  vibration 
of  shells.  This  includes  also  structures  which  undergo  a 
near-uniform  nonlinear  axial /membrane  stress,  due  to  boundary 
constraint.  However,  this  type  of  problem  can  often  be 
solved  adequately  with  conventional  methods. 

o  Boundary  constraints  on  stretched  membranes,  plates,  and 

shells  can  cause  rapidly  varying  local  rotations  and  nonlinear 
strains  at  locations  where  the  boundaries  undergo  sharp  shape 
changes  (corners,  etc.).  Hence,  this  type  of  problem  requires 
the  stability  type  of  element  in  cases  where  accurate  stress 
analysis  within  these  zones  is  desired. 
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o  Prediction  of  instability  behavior  in  general  requires  a 
nonlinear-step  type  of  solution  procedure.  The  effects  of 
nonlinear  effects  such  as  mode  switching,  limit  points,  snap- 
through,  and  buckling  influenced  by  prior  information,  are  not 
usually  amenable  to  eigenvalue  analysis.  The  alternative  of 
asymptotic  instability  analysis  involves  very  difficult  calcu¬ 
lations.  For  either  case  a  competent  nonlinear  step  procedure 
is  required  to  obtain  problem  solutions.  The  case  of  follower 
loads  also  falls  In  this  category. 

o  The  important  problems  of  nonlinear  oscillations  (e.g.,  limit 
cycle  predictions)  are  not  generally  solvable  analytically. 

The  finite  element  approach  with  a  competent  nonlinear  dynamic 
solution  procedure  probably  offers  the  only  practical  approach 
to  this  problem.  This  approach  can  evaluate  nonlinear  respon¬ 
ses  for  the  "almost  periodic"  case  as  well  as  the  true  peri¬ 
odic  case,  and  thus  provide  much  information  about  dynamic 
behavior  and  potential  large  amplitude  dynamic  responses  of 
nonlinear  structures. 

2.2  Technical  Method  Descriptions 

This  section  outlines  briefly  some  of  the  details,  and  prop  'sed  modifica¬ 
tions,  of  the  technical  methodologies  to  be  merged  in  this  contract: 
the  stability  elements;  and  the  static  perturbation  nonlinear  stepwise 
method,  as  applied  to  static  and  dynamic  problems. 

Stability  Elements:  The  present  computer  program  (References  2,  3)  for 
the  stability  elements  (hereinafter  called  HMN  elements,  as  in  these 
references)  has  demonstrated  superior  accuracy,  as  compared  to  conven¬ 
tional  elements,  for  the  case  of  large  bending  deformations.  In  addi¬ 
tion,  the  original  work  of  Haftka,  Mallett,  and  Nachbar  (Reference  1) 
showed  that  a  marked  accuracy  gain  was  obtained  from  the  stability- type 
of  beam  element  in  application  to  buckling  solutions  for  beam-columns. 

The  basic  cause  of  the  accuracy  improvement  gained  from  the  stability 
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and  HMN  elements  is  that,  due  to  the  "strain-smoothing"  enforced  on  the 
nonlinear  strains  by  the  high  order  membrane/axial  displacement  func¬ 
tions,  the  elements'  strain  energies  are  reduced  to  near  minimum  values, 
consistent  with  the  magnitude  of  the  overall  deformation  state.  Since 
the  improvement  is  effected  through  the  membrane  strains,  to  which  cor¬ 
respond  very  large  stiffness  terms,  the  accuracy  gain  can  be  very  large. 
In  the  case  of  stepwise  linear,  nonlinear  problem  solutions,  the  gain  is 
effected  through  the  residual  load  magnitudes.  In  the  case  of  eigen¬ 
value  solutions,  it  occurs  in  the  eigenvalue  itself. 

The  extension  from  one-dimensional  (Reference  1)  to  two-dimensional  ele¬ 
ments  (References  2,  3)  creates  many  difficulties  in  applying  the  origi¬ 
nal  stability  element  concepts.  This  difficulty  resides  primarily  in 

* 

the  fact  that  the  added,  higher  order,  membrane  displacement  functions 
(the  basic  approach  of  the  stability  elements)  are  nonzero  over  the 
entire  two-dimensional  element,  including  its  boundaries.  If  one 
attempts  to  minimize  the  strain  energy  on  the  elemental  level,  which 
would  be  a  relatively  simple  task,  in  general  inter-element  displacement 
incompatibilities  will  bv  created.  The  alternative  is  to  derive  speci¬ 
fic,  explicit  constrains  on  the  added  functions,  such  that  specific 

higher  order  terms  in  the  strains  (e  ,  e  ,  y  )  are  set  to  zero,  without 

x  y  xy 

violating  inter-element  compatibility.  This  alternative  becomes  very 

complicated,  but  nevertheless  was  the  one  adopted  in  the  HMN  element 

work  of  references  2  and  3.  The  work  was  very  successful  for  large 

bending  deformations,  and  less  so  for  large  torsional  deformations.  The 

reason  for  this  is  that  the  specific  higher  order  membrane  functions 

which  compensate  for  large  torsion  (nonlinear  y  }  may  in  some  cases 

xy 

cause  undesirable  higher  order  direct  strains  e  and  e  .  The  reouire- 

x  y 

ment  to  allow  arbitrary  element  orientations  relative  to  any  structural 
deformation  pattern  causes  this  difficulty  to  go  both  ways:  HMN  com¬ 
pensations  for  large  torsion  may  create  undesirable  direct  strains;  HMN 
compensation  for  large  bending  may  cause  undesirable  shear  strains.  The 
physical  meaning  of  this  situation  is  that  either  large  bending  or  large 
torsion  will  in  actual  practice  cause  a  trade  to  occur  between  higher 
order  shear  and  direct  membrane  strains,  such  that  the  structural  poten¬ 
tial  energy  is  minimized.  The  failing  of  the  HMN  elements  of  references 
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2  and  3  is  that  they  deal  with  the  strains  separately,  rather  than  with 
the  total  deformation  state. 

There  are  several  alternatives  for  continuing  work  on  the  present  ele¬ 
ments.  First,  it  is  recognized  that  they  perform  well  as  they  are  cur¬ 
rently  formulated.  They  might  perform  better  with  the  torsion-membrane 
shear  interactions  removed,  which  would  be  very  simple  to  accomplish. 
Finally,  a  method  for  obtaining  the  shear-direct  strain  "trade"  could  be 
devised  and  implemented.  It  appears  that  before  any  of  these  alterna¬ 
tives  are  pursued,  another  option  should  be  investigated.  Figure  1 
shows  an  isoparametric  quadrilateral  (of  the  general  type  of  Reference 
9)  which  has  a  special  relationship  between  nodes  and  displacement 
freedoms.  The  element  has  17  nodes,  of  which  only  8  nodes  are  used  to 
define  the  bending  freedoms,  and  all  17  are  used  to  define  the  membrane 
freedoms.  This  element  will  have  higher  order  membrane  strains,  to 
compensate  the  nonlinear  strains  due  to  bending  and  torsion,  by  virtue 
of  its  extra  9  membrane  only  nodes.  Thus  it  is  basically  a  stability 
element  in  the  sense  defined  by  reference  1.  The  element  has  an  advan- 
age  over  the  HMN  elements  of  references  2  and  3  because  its  higher  order 
freedoms  are  nodally  defined,  and  thus  can  be  committed  to  the  global 
solution  process  without  creating  inter-element  incompatibilities.  The 
displacement  functions  for  the  8  node  bending  behavior  will  be  those  of 
references  2  and  9.  Those  for  the  17  node  membrane  functions  will 
follow  the  conventional  forms  for  isoparametric-elements.  It  is  proposed 
to  use  this  element  in  the  research  described  herein. 

Figure  1  also  shows  a  beam  element  which  will  be  developed.  This 
element  differs  from  conventional  beam  (cubic  displacement)  elements. 

It  has  identical  displacements  to  those  of  one  side  of  the  quadrilat¬ 
eral.  This  will  make  the  two  elements  nodally  compatible  in  problem 
solutions. 

The  work  of  references  2  and  3  includes  many  features  which  are  not 
dependent  on  the  explicit  strain  constraints  of  the  HMN  elements.  These 
include  the  developed  solution  procedure  details,  geometrical  transforma¬ 
tions,  and  nonlinear  shell  equations.  All  of  these  are  applicable  to 
the  element  of  Figure  1,  and  will  be  retained  In  the  proposed  work. 
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The  updating  (or  convection)  of  the  element  coordinate  system  is  pres¬ 
ently  done  for  every  iteration  of  every  solution  step.  This  is  costly, 
and  is  not  always  necessary,  as  this  set  of  transformations  is  only 
important  when  significant  rotations  have  occurred  during  the  step.  It 
is  proposed  to  make  this  updating  conditional  on  the  rotation  magni¬ 
tudes.  The  residuals  will  be  referred  to  the  start-of-step  coordinate 
system  unless  the  updating  is  found  to  be  required.  In  addition,  it 
appears  that  when  it  is  necessary  (rotations  are  large)  to  update  the 
element  coordinate  system,  also  the  solution  coordinate  systems  and  the 
stiffness  matrix  should  be  updated.  The  programs  have  this  feature 
.already  and  it  is  simply  necessary  to  make  the  implementation  condi¬ 
tional  on  the  coordinate  system  updating.  The  changes  to  be  made  will 
cause  the  updating  to  be  done  infrequently,  conditional  on  the  rotation 
magnitudes  being  of  the  order  of  15°.  This  will  reduce  costs  consider¬ 
ably  without  degradation  of  accuracy. 

Several  features  of  the  present  nonlinear  element  formulations  which 
have  proved  particularly  effective  and  will  be  retained  are  listed 
below: 

o  The  iteration  procedure  which  alternates  axial /membrane  and 
all -freedom  iterative  corrections  will  be  retained  (unless  it 
is  shown  to  be  unnecessary  due  to  the  use  of  the  nonlinear 
step  solution  procedure). 

o  The  conditional  updating  of  the  geometric  stiffness  matrix, 
based  on  the  magnitudes  of  the  residual  stresses,  will  be 
retained. 

o  The  convected  coordinate  system  approach  will  be  retained. 

o  The  user-option  of  over-riding  the  internally  computed  solu¬ 
tion  coordinate  systems  is  needed  for  generality  of  boundary 
condition  specification  in  the  nonlinear  case,  and  will  be 
retained. 
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The  shallow-shell  formulation  will  be  retained. 


Further  development  of  the  triangular  HMN  element  is  not  proposed  herein. 
This  element  has  the  recognized  difficulty  of  inter-element  bending 
slope  discontinuities.  While  this  effect  is  not  always  necessarily  a 
bad  one,  it  has  complicated  the  handling  and  interpretation  of  residual 
loads  in  the  stepwise  solution  of  nonlinear  problems.  It  is  noted  how¬ 
ever,  that  the  triangular  element  appears  to  be  nearly  free  of  the 
difficulty  regarding  bending/torsion  and  shear/direct  strain  inter¬ 
actions  which  are  described  above.  Thus,  it  may  ultimately  turn  out 
that  the  triangular  HMN  (BCIZ-Reference  8)  element  merits  further  work. 

Nonlinear  Step  Static  Solution  Procedure:  The  nonlinear  step  capability 
will  be  developed  based  on  the  "Static  Perturbation"  method.  This 
method  was  described  by  Sewell  (Reference  4.)  and  extended  in  a  cost 
effective  manner  to  finite  element  applications  by  Vos  (References  5, 

6).  In  this  procedure  the  nodal  displacement  vector  is  expressed  in 
Taylor  series  form  in  terms  of  a  path  parameter.  Displacement  deriva¬ 
tive  vectors  for  use  in  the  Taylor  series  are  determined  from  solutions 
of  successive  differentiations  of  the  equilibrium  equations,  using  the 
system  tangent  stiffness  matrix.  Problem  solutions  are  determined  from 
the  Taylor  series  expansion.  The  residual  load  method  is  still  used  to 
assure  close  conformance  to  the  equilibrium  path. 

This  nonlinear  step  approach  allows  solutions  to  be  continued  through 
limit  point  instabilities.  The  method  can  incorporate  both  material  and 
geometric  nonlinearities,  as  well  as  the  effects  of  nonconservative 
follower-type  forces.  The  only  matrix  decomposition  required  is  that  of 
the  system  stiffness  matrix,  and  this  is  only  required  once  per  step. 
Techniques  will  be  developed  for  selecting  appropriate  step  sizes.  It 
is  proposed  that  both  second  order  (quadratic  step)  and  third  order 
(cubic  step)  approaches  be  Incorporated  and  compared  for  relative  effi¬ 
ciency.  Appendix  A  gives  the  basic  equations  of  the  static  perturbation 
method  for  J,.z  case  in  which  quadratically  varying  in  step  solution 
variables  are  retained.  Appendix  B  gives  formulas  for  the  nonlinear 
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stiffness  matrices,  in  terms  of  element  displacement  forms  and  material 
property  matrices. 

Figure  2  shows  results  computed  for  a  simple  nonlinear  problem,  comparing 
quadratic  static  perturbation  solutions  and  Newton-Raphson  (piecewise 
linear)  solutions  for  two  step  sizes.  The  static  perturbation  procedure 
Is  seen  to  converge,  with  decreasing  step  size,  much  faster  than  the  con¬ 
ventional  piecewise  linear  method.  Also  shown  is  a  result  computed  with 
the  static  perturbation  method  using  automatically  varied  step  sizes, 
computed  during  the  solution  by  the  formula  aS  =  constant  x  (Q/Q).  The 
results. are  excellent.  The  figure  notes  the  numbers  of  steps  computed 
for  each  plotted  curve. 

For  use  with  the  corrected  coordinate  system  procedure  (updated  total - 
Langrangian  formulation),  the  static  perturbation  method  must  accomplish 
coordinate  transformations  on  the  in-step  nonlinearity  matrix  (^  -  see 
Appendix  A).  The  proposed  method  for  accomplishing  this  is  as  follows: 

Solution  variable  rates  are  computed  in  solution  or  global  system: 

0 

Transform  to  element  baseplane  system 

Olq 

Evaluate  elemental  matrix  PI k*q 

(see  Appendix  A,  Equations  A4,  A5) 

Transform  to  solution  or  global  coordinates 

PI k»q  I  P1K-0 

Form  PI  »  (P1K-0)  0 

This  procedure  avoids  the  requirement  to  transform  the  third  order  tensor 
quantity  PI k.  The  transformation  of  Plk*q  is  a  simple  stiffness  matrix 
transformation,  using  a  conventional  coordinate  transformation  matrix,  T. 
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Nonlinear  Step  Dynamic  Solution  Procedure:  There  are  many  different 
forms  of  discrete  step  solution  procedures  in  use  for  solving  transient 
dynamic  analysis  problems.  For  the  most  part  these  methods  are  based  on 
using  the  set  of  previously  computed  solution  steps,  together  with  the 
differential  equations  of  motion,  to  predict  the  solution  values  at  the 
end  of  the  current  computation  step.  The  stepwise  equations  used  are 
based  on  either  difference  formula  representation  of  timewise  deriva¬ 
tives  of  the  unknown  variables  (using  past  and  future  solution  sets),  or 
on  interpolation  formula  representation  of  these  variables  (again  using 
past  and  future  solution  sets)  with  corresponding  analytical  representa¬ 
tions  of  the  time  derivatives.  In  all  cases  the  equilibrium  equations 
are  forced  to  be  satisfied,  in  terms  of  solution  variables  at  discrete 
time  points,  at  a  particular  point  in  time.  The  choice  of  this  time 
point  is  such  that  the  unknowns  to  be  aetermined,  i.e.,  displacements  at 

c  t 

the  (n+1)—  time  point,  appear  in  the  discretized  equations.  The  dif¬ 
ference  formula  and  interpolation  formula  approaches  are  closely  related, 
but  in  general  lead  to  different  equations,  and  hence  to  somewhat  differ¬ 
ent  numerical  results  in  applications.  Other  distinctions  between  these 
methods  include  whether  the  equations  are  implicit  (solution  requires 
iterations  at  each  time  point,  because  equation  coefficients  are  depend¬ 
ent  on  future  points),  cr  explicit  (solution  steps  do  not  require  itera¬ 
tion  because  equation  coefficients  are  only  dependent  on  past  points); 
and  also  what  order  of  derivatives  are  employed  in  the  equations  of 
motion.  Regarding  the  latter  options,  one  can,  for  example,  simply  use 
the  second  order  equation  of  motion, 

Mflf *  P  -  Ks  Q 

or  employ  further  differentiations  to  obtain,  in  addition. 


In  these  equations,  Kg  and  Ky  are,  respectively,  the  structural  secant 
and  tangent  stiffness  matrices.  A  further  important  consideration 
relative  to  these  solution  procedures  is  whether,  at  each  solution  step 
(and  perhaps  at  each  iteration  of  implicit  methods),  the  computations 
require  only  the  decomposition  of  the  structural  mass  matrix,  or  alter¬ 
natively,  a  decomposition  involving  the  mass  and  stiffness  matrices. 

The  latter  is  generally  the  case  for  implicit  methods,  and  is  very 
costly  in  practical  numerical  work. 

The  implicit  methods  in  some  instances  have  the  advantage  of  uncondi¬ 
tional  stability  as  the  time  steps  are  increased  in  size, ‘while  the 
explicit  methods  become  unstable  for  particular  step  sizes  (of  the  order 
of  the  half-period  of  the  highest  frequency  components  of  the  structural 
system).  The  advantage  of  the  unconditional  stability  is  that  the 
highest  frequency  structural  actions  of  a  finite  element  model  (which 
can  be  of  very  high  frequency  for  fine  discretizations)  will  be  "damped" 
to  a  near  zero  amplitude  in  problem  solutions.  However,  particularly 
for  nonlinear  problems,  obtaining  good  solution  accuracy  may  require 
smaller  time  steps  for  properly  representing  rapidly  varying  structural 
behavior  than  would  be  required  to  satisfy  stability  criteria  for  the 
integration  procedure.  Thus  it.  appears  that  the  implicit  methods, 
requiring  costly  stiffness  matrix  decomposition,  may  not  be  optimum  for 
nonlinear  dynamic  analysis.  In  addition,  the  implicit  methods  impart  a 
numerically-induced  artificial  damping  to  problem  solutions,  which  in 
itself  requires  the  use  of  small  time  steps  to  avoid  excessive  energy 
loss  due  to  the  artificial  damping  effects. 

References  7,  12-16  discuss  various  solution  procedures  of  the  general 
types  described  above.  The  discussions  In  these  references  are  for  the 
most  part  mathematical  in  approach.  In  order  to  put  such  methods  in 
perspective,  a  particular  procedure,  called  the  Houbolt  method  (Reference 
16)  and  generally  considered  to  be  a  superior  method,  has  been  used  to 
solve  a  simple  nonlinear  problem.  Figure  3  shows  the  numerical  results 
for  several  time  step  sizes.  The  solution  involves  Iterations  at  each 
time  point,  and  the  data  shown  are  iterated  to  obtain  fully  converged 
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results.  It  will  be  seen  in  later  discussions  that  the  accuracy  of  the 
Houbolt  method,  at  least  in  this  particular  nonlinear  problem,  is  not 
particularly  good,  and  can  be  improved  on  by  simpler  methods.  The 
Houbolt  method  is  described  in  Appendix  A. 

A  basically  different  type  of  formulation  starts  from  the  representation 
of  the  solution  as  a  Taylor  series.  In  this  case  the  solution  at  the 
(n+1)—  time  point  is  based  on  its  derivatives  at  the  nth  time  point. 

This  approach  offers  a  number  of  advantages:  complete  freedom  to  vary 
time  step  size  during  the  solution;  solution  behavior  governed  by  the 
most  recent  structural  behavior,  rather  than  by  extrapolation  from  past 
behavior;  simple  extension  to  higher  orders  of  approximation,  even  dur¬ 
ing  a  problem  solution,  without  changing  the  basic  solution  equations; 
ability  to  handle  in-step  nonlinearity  without  the  use  of  an  implicit/ 
iterative  solution  method  (only  the  mass  matrix  needs  to  be  decomposed).  . 
This  approach  is  analogous  to  the  static  perturbation  procedure,  and  the 
relevant  equations  are  given  in  Appendix  A.  This  approach  is  proposed 
for  the  subject  research  and  computer  program  development. 

The  Taylor  series  representation  approach,  called  herein  the  "dynamic 
perturbation  method",  can  be  formulated  to  make  use  of  the  second  order 
equilibrium  equation,  plus  an  arbitrary  succession  of  higher  order  equa¬ 
tions  obtained  by  differentiating  the  basic  equation.  Through  the 
higher  order  derivatives,  more  complete  information  describing  the 
variation  of  the  forces  acting  during  the  computation  time  step  is 
incorporated  into  the  solution.  This  is  clearly  seen  in  Equations  1-3, 
in  which  KjO  represents  the  effect  of  variable  force  at  constant  stiff¬ 
ness,  and  fry  0,  is  the  first  term  which  represents  the  effect  of  in-step 
structural  nonlinearity.  Equations  1-3  can  be  solved  for  Q,  Q,  Qiv,  •••• 
etc.,  requiring  only  decomposition  of  the  mass  matrix,  M.  These  deriva¬ 
tives,  evaluated  at  time  tn  are  used  in  the  Taylor  series  (about  tp) 
through  which  the  solution  at  time  tp+1  is  computed.  The  simplest 

option,  using  only  Q,  does  not  generally  provide  accurate  problem  solu- 
•••  •••  H  v 

tions.  Including  Q,  or  Q  and  Q  ,  causes  the  results  to  be  very  accurate, 
even  for  time  steps  approaching  the  stability  limit  {ht\  1/2  period)  of 
the  formulation. 
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The  Inadequacy  of  the  lowest  order  "dynamic  perturbation"  method  is 
easily  remedied  by  a  slight  change  in  formulation,  described  in  Reference 
15.  The  resulting  particularly  simple  method,  called  the  "acceleration- 
pulse"  method  (Reference  15),  offers  probably  the  most  cost  effective  of 
the  available  solution  procedures.  This  method  achieves  its  excellent 
accuracy  by  compensating  errors,  as  described  in  Appendix  A.  Since  it 
is  only  a  zero'th  order  method  (based  on  Q  only),  only  simple  calcula¬ 
tions  are  involved,  and  the  method  does  not  include  effects  of  in-step 
nonlinearities.  Nevertheless,  because  of  its  ease  of  use,  economy  and 
good  accuracy,  it  is  felt  that  this  method  should  be  included  in  the 
subject  program  development,  and  it  is  proposed  to  include  it  as  a  user- 
optional  choice,  along  with  the  Taylor  series,  or  "dynamic  perturbation", 
method. 

Figure  4  illustrates  the  acceleration  pulse  method  and  the  "dynamic 
perturbation"  method  through  the  4th  derivative  for  a  simple  problem. 

In  this  problem,  the  approach  keeping  w^v  is  essentially  exact,  as 
proved  by  solutions  obtained  with  a  set  of  smaller  time  steps.  The  data 
illustrate  elastic,  plastic,  material  failure,  and  load  discontinuity 
induced  behavior.  The  superiority  of  the  higher  order  method,  which 
includes  both  linear  and  nonlinear  in-step  force  variation,  is  seen  to 
be  greatest  when  some  degree  of  discontinuity  of  load  or  stiffness 
behavior  is  present,  particularly  when  the  discontinuity  is  an  added, 
positive  load.  Even  in  this  cas,,  however,  the  excellent  accuracy  of 
the  acceleration  pulse  method,  in  relation  to  its  simplicity,  is  clearly 
seen.  The  Houbolt  method  (Figure  3}  was  seen  to  provide  mediocre 
results  in  comparison  with  the  "dynamic  perturbation"  method,  even  for 
the  simple  elastic  case. 

It  should  be  noted  that  the  simple,  one-degree-of-freedom  example  may  be 
somewhat  misleading.  Judgement  suggests  that  more  severe  calculation 
difficulty,  with  attendant  greater  accuracy  requirements,  would  be  pres¬ 
ent  in  multi -degree-of-freedom  problems,  particularly  when  material 
yield  or  failure  occurs,  resulting  in  growth  and  contraction  (unloading) 
of  failure/yield  zones.  Difficulties  related  to  this  type  of  behavior 


Cl  4 


have  previously  been  encountered  with  the  "acceleration-pulse"  method 
(Reference  17). 

Reference  7  uses  an  implicit,  interpolation-type,  solution  procedure 
retaining  u*  to  solve  Inelastic  problems  of  beams  and  shells.  The  method 
is  related  to  the  Houbolt  method.  Despite  the  complexity  and  inherent 
cost  of  the  method,  small  time  steps  were  apparently  required  to  obtain 
accurate  solutions.  This  may  suggest  that  some  sort  of  "dynamic  residual 
load"  concept  would  be  a  valuable  asset  with  this,  and  probably  other, 
solution  procedures.  Such  a  residual  load  method  will  be  investigated 
as  an  option  In  the  proposed  computer  program. 

The  goal  of  the  proposed  research  is  to  handle  nonlinear  dynamic  problems 
with  relatively  large  time  steps  (of  the  order  permissible  for  linear 
problems,  governed  by  solution  stability  criteria),  while  using  a  solu¬ 
tion  procedure  which  only  requires  decomposing  the  mass  matrix.  The 
latter  assures  a  method  which  Is  both  fast  and  simple.  The  approaches 
proposed  (dynamic  application  of  static  perturbation  procedure,  and  the 
acceleration  pulse  method,  Appendix  A)  provide  these  desirable  features. 
In  addition,  the  first  method  lends  itself  to  the  automatic  computation 
of  time  step  size,  based  on  specified  accuracy  criteria  (using  ratios  of 
time  derivatives  of  solution  quantities).  In  most  problems,  this  can 
yield  considerable  savings  in  computing  costs. 
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APPENDIX  A 


This  appendix  briefly  describes  the  equations  of  the  numerical  stepping 
procedures  (static  and  dynamic)  considered  in  this  proposal.  In  the 
equations  given,  the  following  definitions  hold  (matrix  notation  omitted). 

Q  =  solution  vector  (col urn  matrix) 

P  =  load  vector  (column  matrix) 

M  *  mass  (square  matrix) 

=  secant  and  tangent  stiffness  (square  matrix) 

C  =  damping  coefficient  (square  matrix) 

At  =  incremental  time  or  incremental  path  parameter 

PI  *  load  vector  which  accounts  for  in-step  internal  structural 
loads  due  to  nonlinearity  (column  matrix) 

PI K  =  the  rate  of  change  of  Ky,  due  to  nonlinearity  (third  order 
tensor,  or  "cubic  matrix  array") 

(’),(  )»etc.  *  denotes  time  or  "path  parameter"  derivatives 

(  )  *  denotes  the  nth  time  point  or  path  parameter  point 

Static  Perturbation  Method  (through  Q) 

The  starting  equation  is  the  equilibrium  equation 

P  =  KSQ  A1 
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Successive  differentiations  yield 


P  -  K$Q  +  K$Q  =  KtQ  A2 

M  **  •  • 

P  =  ICj-Q  +  KtQ  =  ICj-Q  +  PI  A3 

where  PI  can  be  written  as 

•  • 

PI  =  (P1K-Q)  Q  A4 

•  • 

and  Kt  =  PTK-Q  A5 

Solving  for  the  derivatives  of  Q, 

Q  s  Kj1  P  A6 

•.  * 

Q  *  (*P  -  PI )  A7 


It  is  noted  that  only  needs  to  be  inverted  (decomposed),  even 
though  the  equations  contain  the  effects  (through  PI)  of  structural 
nonlinearity. 

The  final  solution  is  obtained  by  a  Taylor  series  stepping  process 
in  which  Qn+1  is  computed  from  the  previous  step  solution  Qn  and  the 
start-of-step  derivatives  (Equations  A6,  A7)  Qn ,  Qn 

<U  ■  Q„  +  %  4t  *  1/2  (lt)2  AB 
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The  procedure  can  be  used  retaining  only  Qn,  in  which  case  it  is  equiva¬ 
lent  to  the  simple,  and  usually  inadequate,  piecewise  linear  Newton- 
Raphson  procedure.  The  real  accomplishment  of  the  static  perturbation 

methods  lies  in  including  the  higher  derivatives.  It  is  noted  that,  by 
*  «  •  • 

retaining  PI,  Q  can  be  included,  and  similarly  even  higher  derivatives 
can  easily  be  included.  See  Appendix  B  for  closed  form  equations  for 
nonlinear  stiffness  matrices. 

Dynamic  Perturbation  Method  (through  Q1v) 

The  starting  equation  is  the  second  order  equation  of  motion,  with  time 
the  path  parameter, 

MQ  =  P  -  K$Q  -  CQ  A9 

Differentiating,  and  solving  for  successive  derivatives. 


Q  *  M"1  £P  -  K$q  -  CQ] 

Q'  =  M-1  [P  -  ICpQ  -  CQ  -  CQ] 


Q1V  =  M-1  [P  -  ICj-Q  -  PI  -  CQ  -  2CQ  -  CQ] 


The  Taylor  series  about  time  tn  gives  the  solution  at  tp+1 

Vi  *  «n  +  %  Lt  +  «n  ^  +  \  ^  +  C  ^  A13 

Vi  -  i  *  %  « +  \  ^  ^  A14 

It  is  noted  that  in  the  dynamic  case,  both  Qn+^  and  Qn+^  are  solved  for, 
in  order  that  the  succeeding  steps  can  be  handled  as  an  initial  value 
problem. 

If  terms  are  only  retained  through  Q,  the  method  is  not. accurate.  The 
physical  reason  for  this  is  that  the  computed  value  of  Qn+1  consider: 


only  a  constant  acceleration  through  the  step  -  i.e.,  a  constant-force 
step.  This  is  not  adequate  for  even  linear  analysis.  If  terms  are 
retained  through  Q,  the  method  has  effectively  retained  in-step  linear 
force  variation,  through  the  term  KjQ.  This  level  of  approximation  has 
been  found  to  be  very  accurate  for  moderately  nonlinear  behavior. 
Retention  of  Q1V  includes  in-step  nonlinearity  and  further  improves 
accuracy  for  strongly  nonlinear  behavior. 

"Acceleration  Pulse"  Method 


This  method  can  be  derived  from  Equation  A9  (without  the  damping  term)  by 
using  a  central  difference  formula  for  Q,  and  representing  the  start-of- 
step  velocity  Qn  by  a  backward  difference  formula.  The  result  is  equiva¬ 
lent  for  a  rather  surprising  modification  of  Equations  A13  and  A14,  as 
follows: 


q;+i  *  «Vi  -  «»>'«  A16 

The  starred  quantities  indicate  approximate  velocities.  The  appearance 
of  the  extra  acceleration  term  in  equation  A15  compensates  for  the  error 
incurred  by  the  backward  difference  representation  of  the  velocity  in 
Equation  A16.  It  can  be  rather  easily  seen  that 

Qn  *  Q*  +  5- Qn(^)2  A17 

with  the  result 

Vl  ■  +  «  +  7  <«>’  5n  A18 


A19 
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Equation  A18  differs  from  Equation  At 3  retaining  only  through  Q  in  that 
Qn  in  A18,  obtained  from  Equation  AT 7,  is  much  more  accurate  than  the 
corresponding  term  in  Equation  A13. 

The  acceleration  pulse  method  achieves  truly  outstanding  accuracy,  in 
consideration  of  its  simplicity,  even  for  quite  large  step  sizes.  The 
means  of  including  damping  while  maintaining  the  internal  error  compensa¬ 
tion  feature  of  the  method  is  not  developed  as  yet. 


Houbolt  Method 


joint  backward  difference) 


The  Houbolt  method  uses  a  four  point  interpolation  formula  for  Q,  based 
on  the  unique  cubic  polynomial  passed  through  four  equally  spaced  points 

5  *  V3  +  St  T  V3  +  3C|ti-2  ‘  I  Vl  +  T  V 

+  (  St  >2  tQn-3  '  I  +  2Vl  +  ?  V 
*  <st  >’  fl2° 

Differentiating  this  formula  to  obtain  Qn+,  and  Qn+1 .  and  substituting  in 
Equation  A9  yields 

[M  +  "  6t.c  ♦  1.  (it )2, K$  ]  Qn+1  .  J.  !«)!Pn+1 
n+1  n+1 

+  Qn  (  |m  +  |  A C-C)  -  Qn-1  (2M  +  |  At-C) 

+  Qn_2  (  \  M  +  ^  At-C)  A21 

Use  of  this  equation  requires  inversion  (decomposition)  of  the  quantity 

[M  +  j^-At  Cn+1  +l(At)2Ks  ], 

n+1 


which  makes  the  method  implicit,  as  K$  is  not  generally  known  until 
Qn+1  ’s  *(nown*  Solutions  are  obtainednS}  iteration. 
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APPENDIX  B 


This  appendix  briefly  describes  the  equations  required  for  the  elemental 
level  forces  and  stiffness  matrices,  and  their  derivatives,  for  nonlinear 
problems.  In  the  equations  given,  the  following  definitions  hold  (matrix 
notation  omitted): 

q  =  element  freedoms  (column  matrix) 

a  =  stresses  or  stress  resultants  (column  matrix) 

e  »  Lagrangian  strains  or  curvatures  (column  matrix) 

9  =  spatial  derivatives  of  displacements  (column  matrix) 

D  »  material  stress/strain  relation  (square  matrix) 

AO  =  linear  terms  for  Lagrangian  strain  definitions  (rectangular 
matrix) 

A1  *  nonlinear  terms  for  Lagrangian  strain  definition  (3rd  order 
tensor,  or  "cubic  matrix  array") 

A  =  terms  for  nonlinear  Lagrangian  strain  rate  definition 
(rectangular  matrix) 

B  *  strain/displacement  rate  relation  (rectangular  matrix  ) 

G  *  element  shape  function  derivatives  (rectangular  matrix) 
kT  =  element  tangent  stiffness  (square  matrix) 
pi  =  element  nonlinear  structural  force  terms  (column  matrix) 

(•).  ("),  etc.,  denote  time  or  "path  parameter"  derivatives 
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The  equations  for  elemental  level  force  and  stiffness  quantities  can  be 
derived  in  a  straightforward  manner,  based  on  virtual  work.  The  basic 
relations  are  provided  here.  The  stress/strain  relation  is  given  by 

o  «  D  e  B1 

This  equation  is  formulated  so  as  to  include  treatment  of  sandwich  and 
nonisotropic  materials.  The  quantities  in  B1 ,  as  well  as  the  element 
displacement  (shape)  functions  and  displacement  derivatives,  are  evaluated 
at  a  series  of  numerical  integration  points  within  the  element.  The 
displacement  derivatives  9  are  given  by 

6  *  G  q  B2 

Strains  are  defined  by 

e  a  (AO  +  1/2  A1  6)  0  B3 

and  strain  rates  by  * 

e  =  (AO  +  A1  8)  §  =  A  0  B4 

Combining  B2  and  B4  provides 

e  =  A  G  q  =  Bq  B5 

The  virtual  work  formulation  leads  to  the  expression  for  element  nodal 
forces,  in  terms  of  a  numerically  integrated  volume  integral 

p  =  G  A  o  dV  =  Jv  B  a  dV  B6 

Substituting  B1  and  B3  into  86,  and  differentiating,  provides  the  first 
order  (tangent  stiffness)  relation 


p  ■  kT  q 


B7 


where 


kT  *  ;v  G(oAl  +  ADA)  G  dV 


B8 


Here  kj  is  a  symmetric  matrix  due  to  symmetry  properties  of  A1 .  Differen¬ 
tiating  B7  provides  the  second  order  relation  for  the  element 


p  *  kT  q  +  q  =  kT  q  +  pi 


B9 


where 


pi  *  kT  q 


BlOa 


In  an  expanded  but  perhaps  more  computationally  advantageous  form 


pi  *  'v  GD  (2AU0  +  AAl86)dV 


B10B 


d 
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